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1.0.0.  Introduction 

The  basic  problem  which  Is  being  addressed  here  Is  that  of  constructing 
a  smooth  (at  least  continuous  first  partial  derivatives)  bivariate  function, 

F(x,  y),  which  takes  on  certain  prescribed  values,  F(xk,  yk)  -  fk,  k  *  1,  ...,  N. 
The  points  (xk,  yk)  are  not  assumed  to  satisfy  any  particular  conditions  as 
to  spacing  or  density,  hence  the  term  "scattered".  It  is  usually  convenient 
to  think  of  the  values  fk  as  arising  from  some  underlying  (not  necessarily 
known)  function  f(x,  y),  so  that  fk  *  f(xk,  yk),  k  ■  1,  ....  N. 

The  problem  of  Interpolation  of  scattered  data  In  two  or  more  Independent 
variables  has  been  addressed  by  numerous  authors,  as  can  be  seen  by  the 
bibliography.  Many  of  the  basic  Ideas  Involved  are  discussed  in  two  survey 
papers  (both  over  a  wider  class  than  we  consider  here)  due  to  Schumaker  [49] 
and  Barnhill  [4].  Some  of  the  Ideas  seem  to  be  mainly  that,  Ideas,  with  only 
a  few  numerical  examples  given,  often  not  well  thought  out  or  very  definitive 
In  terms  of  the  actual  capabilities  of  the  method.  In  addition,  most  of  the 
methods  Involve  one  or  more  ad  hoc  assumptions  requiring  a  user  to  specify 
parameters  (one  or  more).  Generally  only  cursory  attention  has  been  paid  to 
appropriate  choice  of  these  parameters  and  their  overall  effect  on  the  Inter¬ 
pol  ant  has  usually  not  been  determined. 

Out  of  this  situation  arose  a  desire  to  attempt  to  answer  a  number  of 
questions,  basically  all  related  to  the  question:  Which  of  these  many 
methods  deserve  further  study  and  development,  and  which  should  be  discarded? 
Included  here  is  the  determination  of  some  default  values  for  ad  hoc  parameters 
In  methods  which  require  them.  The  default  values  should  give  reasonably  good 
results  over  a  number  of  different  sets  of  data,  and  preferably  the  interpolant 
should  be  rather  stable  with  respect  to  changes  In  the  parameter.  Additionally, 


It  Is  convenient  for  the  user  If  the  parameter  Is  related  to  something 
about  the  data  which  can  be  easily  estimated.  In  many  cases  (perhaps  all), 
subjective  judgements  must  be  made  about  these  matters,  although  some  hard 
Information  can  be  obtained. 

Some  previous  fairly  estensive  work  had  been  done  by  McLain  [39]  which 
Inspired  a  Somewhat  similar  study  of  another  class  of  Ideas  by  the  current 
Investigator  [16].  The  Initial  thrust  of  the  Investigation  was. to  compare  a 
few  "local"  methods  to  determine  which  seem  to  work  reasonably' well .  As 
the  Investigation  proceeded,  more  Ideas  were  supplied  by  colleagues  and 
others,  so  that  In  the  end,  more  than  a  few  methods  are  tested  and  compared 
here.  Including  "global"  methods.  The  total  number  of  programs  involved 
In  this  study  Is  29,  some  of  which  are  fairly  minor  variations  of  others. 

The  concept  of  a  "global"  method  Is  easily  understood.  The  Interpolant 
Is  dependent  on  all  data  points,  and  addition  or  deletion  of  a  data  point, 
or  a  change  of  one  of  the  coordinates  of  a  data  point  will  propagate  through¬ 
out  the  domain  of  definition.  The  Idea  of  a  "local"  method  Is  not  so  clear. 
Typically  one  thinks  of  It  meaning  that  addition  or  deletion  of  a  point,  or 
a  change  of  one  of  the  coordinates  of  a  data  point  will  affect  the  interpolant 
only  at  nearby  points,  that  Is,  the  Interpolant  will  be  unchanged  at  distances 
greater  than  some  given  distance.  There  are  some  difficulties  here.  If  the 
data  (the  (x^,  y^)  points)  are  "random",  one  must  Inspect  (In  someway)  all 
the  data  to  determine  which  are  "nearby".  Does  this  mean  there  Is  no  such 
thing  as  a  "local"  method?  (Rosemary  Chang  first  mentioned  this  Idea).  We 
have  taken  a  somewhat  more  liberal  view  of  "local"  and  take  It  to  mean  that 
the  Interpolant  Involves  only  "nearby"  points  and  one  or  more  parameters. 

We  allow  the  parameters  to  have  been  globally  determined  as  a  matter  of  user 
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convenience,  even  though  a  (successful)  argument  can  be  made  that  then  the 
method  Is  not  local.  Thus,  we  classify  methods  as  local  or  global  without 
regard  to  how  parameters  are  chosen  or  computed. 

The  use  of  global  methods  Is  not  feasible  for  very  large  N  since  they 
often  Involve  the  solution  of  a  system  of  0(N)  equations  (often  exactly  N) 
and  In  any  case  Involve  processing  all  points.  When  systems  of  equations 
must  be  solved,  the  systems  are  often  full  and  not  well  conditioned.  While 
our  primary  aim  was  to  Investigate  local  methods  suitable  for  very  large  data 
sets  (several  hundred  points  up  to  some  millions,  say),  In  many  Instances  local 
methods  Involve  the  use  of  global  methods  on  smaller  sets  which  are  then 
"blended"  together  to  obtain  a  locally  defined  global  Interpolant.  Thus  it 
makes  sense  to  test  global  methods  on  moderately  sized  sets  of  data.  By 
the  same  token.  It  Is  not  necessary  to  test  local  methods  on  sets  of  10000 
points  (say)  by  virtue  of  the  fact  that  they  are  local.  If  very  large  sets 
of  data  were  to  be  considered,  It  is  clear  that  a  different  Implementation 
approach  might  be  necessary,  one  which  would  involve  a  larger  amount  of  pre¬ 
processing  and  perhaps  additional  storage. 


1.1.0.  Tested  Characteristics  of  Methods 

The  characteristics  on  which  various  methods  are  to  be  compared,  and 
how  they  are  to  be  weighted  In  the  final  analysis,  are  somewhat  subjective 
While  no  representation  Is  made  that  the  list  Is  exhaustive  (or  even  close 
to  It),  nor  that  everyone  will  be  In  agreement  on  It,  the  following  Items 
are  the  ones  considered  here.  We  give  them  and  discuss  them  In  order  of 
decreasing  Importance.  In  the  presentation  of  Information  In  the  summary 
(tables  and  perspective  plots)  each  reader  may  weight  various  aspects  to 
suit  his  own  needs. 
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1.1.1.  Accuracy 

Accuracy  In  reproducing  a  known  surface  Is  certainly  one  Important 
aspect  of  comparison.  In  the  usual  application  no  representation  of  the 
underlying  surface  z  a  f(x,  y)  Is  known,  however.  If  the  method  approximates 
a  variety  of  surface  behavior  faithfully  we  can  expect  It  to  give  reasonable 
results  In  other  Instances.  Quantitative  numbers  can  be  put  on  the  perfor¬ 
mance  of  a  method  tested  In  this  fashion,  and  we  have  used  this  Idea  extensively. 
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1.1.2.  Visual  Aspects 

It  has  developed  during  the  course  of  this  project  that  the  appearance 
of  the  Interpolant  when  viewed  In  a  perspective  plot  Is  very  Important. 
Visual  ratings  are  often  closely  related  to  the  accuracy  with  which  an  Inter 
pol ant  reproduces  test  surfaces,  there  seems  to  be  a  closer  relationship 
when  accuracy  Is  high  since  there  Is  less  chance  for  the  Interpolant  to 
misbehave.  At  moderate  accuracies  one  Interpolant  may  be  visually  pleasing 
while  another  with  similar  accuracy  Is  not. 

The  visual  aspect  Is  quite  subjective  and  ratings  by  different  persons 
will  give  somewhat  different  results,  although  probably  not  contradictory 
ones.  While  It  Is  felt  that  the  visual  aspect  Is  quite  Important,  exactly 
how  this  Information  Is  Integrated  Into  the  overall  assessment  of  a  method 
Is  also  a  subjective  matter,  however  It  Is  rare  that  a  dilemma  occurs  In 
this  study. 
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1.1.3.  Sensitivity  to  Parameters 

Many  of  the  tested  methods  Involve  the  choice  of  one  or  more  parameters. 
These  choices  have  generally  been  converted  to  ones  which  are  related  to 
mean  distances  to  nearest  neighbor,  although  precisely  that  Idea  Is  never 
directly  used.  Here  we  are  talking  of  nearest  neighbor  In  the  set  of  points 
{(xk,  yk)}»  Sometimes  the  parameter  takes  the  form  of  an  anticipated  number 
of  points  In  the  region  which  defines  a  local  Interpol  ant. 

Methods  which  involve  parameters  underwent  informal  testing  for  suitable 
values  of  the  parameters.  Methods  which  survived  this  and  other  tests  have 
parameter  variation  tests  tabulated  In  the  results.  Some  methods  were  found 
to  be  capable  of  generating  creditable  results  for  an  appropriate  value  of 
the  parameter,  but  were  sensitive  to  It,  or  gave  poor  results  on  similar  data 
when  the  same  value  was  used.  These  results  are  mentioned  In  Section  3.  It  is 
desirable  to  have  a  method  which  Is  stable  with  respect  to  changes  In  the 
parameter,  and  such  methods  were  found,  as  we  note  later. 


1.1.4.  Timing 

1  he  computational  effort  required  Is  generally  not  of  great  Interest, 
unless  it  Is  very  high.  In  these  respects,  only  one  method  was  tested  which 
was  discounted  for  this  reason.  Some  methods  are  quite  efficient  in  terms 
of  time  required  for  the  calculations.  These  methods  have  generally  been 
found  deficient  in  other  categories,  unfortunately.  For  methods  which  involve 
a  preprocessing  phase,  distinct  from  an  evaluation  (of  the  interpolant)  phase, 
the  two  times  for  standard  problems  are  given  separately.  Execution  times 
were  taken  from  the  multiprogramming  environment  on  the  IBM  360/67  and  as 
such  may  vary  considerably  with  exactly  the  same  data.  More  is  said  of  this 
later. 
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1.1.5.  Storage  Requirements 

As  with  computational  effort,  storage  requirements  are  not  crucial,  unless 
they  are  very  high.  For  very  large  problems  this  may  be  altered,  of  course. 

We  count  storage  requirements  only  In  terms  of  additional  arrays  needed  to 
store  data  beyond  the  (xk,  yk,  fk)  points.  No  account  is  taken  of  simple 
variables  or  program  length. 
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1.1.6.  Ease  of  Implementation 


Ease  of  •implementation  Is  of  no  great  concern  if  one  obtains  a  working 
program.  In  other  Instances  It  may  be  of  considerable  Importance.  The 
judgement  is  again  subjective.  Further,  it  could  be  different  depending  on 
the  philosophy  behind  the  Implementation.  The  form  of  the  implementation 
could  involve  trade-offs  between  timing  and  storage  and  would  doubtlessly 
alter  the  ease  of  Implementation. 

Implementation  of  programs  specifically  for  this  project  generally  was 
done  with  a  lack  of  frills.  Reasonable  care  was  taken  to  assure  that  a 
grossly  inefficient  algorithm  was  not  coded,  but  no  doubt  It  Is  possible  to 
Improve  on  most  of  them.  In  particular,  use  of  some  preprocessing  and 
additional  storage  was  not  used  to  Increase  efficiency  during  the  evaluation 
phase.  For  a  general  purpose  program  this  should  probably  be  done,  In  many 
Instances.  Some  of  the  documented  programs  did  use  these  devices.  Ease  of 
Implementation  Is  generally  meant  to  take  into  account  the  complexity  of  the 
ideas  Involved  in  the  method  and  the  amount  of  code  required. 
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1.2.0.  The  Testing  Process 

The  Initial  tests  performed  on  a  few  methods  eventually  gave  rise  to  a 
standard  set  of  test  problems  and  a  set  of  supporting  subprograms  to  generate 
statistics  from  the  tests  and  generate  perspective  plots  of  surfaces.  Due  to 
the  evolution  of  Ideas  as  the  study  progressed,  some  aspects  of  the  process 
are  not  as  simple  as  they  might  have  been.  This  Is  particularly  true  of  some 
of  the  test  functions,  but  this  has  no  bearing  on  the  validity  of  the  tests. 


1.2.1.  The  Test  Program 

To  enable  testing  many  different  methods  In  a  consistent  manner,  and  with 
a  minimum  of  effort,  a  set  of  standard  subprograms  was  developed  which 
generate  the  test  cases,  compute  deviation  statistics  for  known  test  surfaces, 
obtain  timing  statistics,  and  generate  and  label  perspective  plots  of  the 
surfaces.  With  the  current  set  of  supporting  subprograms  It  is  generally 
quite  easy  to  test  a  new  method  which  Is  typically  supplied  as  a  subprogram 
(or  several)  which  generates  the  values  of  the  Interpolant  at  a  grid  of  x-y 
points.  Typically  all  that  Is  required  Is  to  set  certain  parameters,  reserve 
any  required  workspace,  and  call  the  subroutine,  all  of  which  can  be  done  with  a 
few  statements  added  to  the  prototype  driver  program. 
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1.2.2.  The  Test  Problems 


The  basic  set  of  test  problems  consisted  of  six  different  test  functions 
over  three  different  x-y  point  sets,  and  two  x-y-z  point  sets  from  the 
literature,  one  of  those  used  In  a  second  version  with  one  of  the  coordinates 
scaled.  Another  Interesting  test  was  the  computation  of  a  "cardinal" 
function  obtained  by  setting  all  function  values  on  a  point  set  to  zero,  save 
one. 

The  six  test  functions  were  all  to  be  approximated  on  [0,  1]  .  Four  of 
them  were  basically  obtained  from  McLain's  paper  [39],  but  were  translated 
to  [0,  1]  from  [1,  10]  and  some  modified  slightly  to  enhance  the  visual 
aspects  of  the  surface.  The  other  two  were  generated  by  the  author  to  provide 
a  fundamentally  different  shape  In  one  case  (saddle),  and  to  provide  a  surface 
with  a  variety  of  behavior  on  one  surface  to  serve  as  a  principal  test  func¬ 
tion. 

The  principal  test  function  Is  given  by 
f-j (x,  y)  -  .75  exp[  -  +  .75  exp[  -  -  ^] 

+  .5  exp[  -  -l—7)2  .  .2  exp[  -  (9x-4)2  -  (9y-7)2]. 

This  surface  consists  of  two  Gaussian  peaks  and  a  sharper  Gaussian  dip 
superimposed  on  a  surface  sloping  toward  the  first  quadrant.  The  latter 
was  Included  mainly  to  enhance  the  visual  aspects  of  the  surface,  which  Is 
shown  In  Figure  4.0.1 .0. 

The  second  test  function,  essentially  obtained  from  McLain  Is 
f2(*,  y)  a  g{tanh(9y  -  9x)  +  1], 

2 

This  surface  consists  of  two  nearly  flat  regions  of  height  0  and  g-,  joined 
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by  a  sharp  rise,  almost  a  cliff,  running  diagonally  from  (0,  0)  to  (1,  1). 
The  test  surface  is  shown  in  Figure  4. 0.2.0. 


The  third  test  function  was  generated  by  the  investigator  and  is 


f,(x,  y) 


1.25  +  cos(5.4y) 
6[1  +  (3x  -  lr] 


This  surface  is  saddle  shaped  and  is  shown  in  Figure  4.0. 3.0. 

The  fourth  test  function,  essentially  obtained  from  McLain,  is 

f4(x,  y)  ■  j  exp[  -  +  (y-|-)  )]. 

This  surface  is  a  Gaussian  hill  which  slopes  off  in  rather  gentle  fashion  in 

o 

[0,  1]  .  It  can  be  seen  in  Figure  4. 0.4.0. 

The  fifth  test  function  was  also  essentially  obtained  from  McLain  and  is 

f5(x.  y)  B  }  exp[  -  JL((x-jr)  +  (y-j)  )]. 

This  surface  is  a  steep  Gaussian  hill  which  becomes  almost  zero  at  the  bound¬ 
aries  of  the  unit  square.  It  can  be  seen  in  Figure  4. 0.5.0. 

The  sixth  test  function  Is  also  essentially  from  McLain,  and  Is 

2  g  — 

f6(x,  y)  ■  ^[64  -  81  ((x-J-)  +  (y-J)  )]2  - 

This  surface  represents  the  part  of  a  sphere  above  the  unit  square.  The 

8  111 

sphere  Is  of  radius  |-w1th  center  at  (£,  -  £) .  The  surface  is  shown  in 

Figure  4. 0.6.0. 

There  were  three  different  sets  of  points  over  [0,  1]  used  in  the  tests, 
The  first  set  consisted  of  100  points  generated  by  a  pseudorandom  number 
generator,  one  point  In  each  square  of  side  ^  centered  at  (^,  jj-)  for 
1,  j  ■  1,  ....  10.  This  yields  a  set  of  scattered  points  forced  to  have 


somewhat  uniform  density,  although  as  can  be  seen  in  Figure  0.1. 0.0.  there 
are  locally  large  variations  In  density.  The  triangulated  set  of  points  Is 
also  shown  In  Figure  0.1. 0.0.  Part  of  the  unit  square  Is  outside  of  the 
convex  hull.  The  points  are  listed  In  Table  1. 

The  second  set  of  data  consists  of  33  points  and  was  generated  by  the 
Investigator  to  purposely  have  some  areas  sparsely  populated  by  points 
while  other  areas  are  not.  This  set  of  points  Is  shown  In  Figure  0.2. 0.0. 

The  points  are  listed  In  Table  2. 

The  third  set  of  points  was  digitized  by  Gregory  M.  Nielson  and  Is 
similar  In  disposition  to  a  set  of  points  appearing  in  McLain  [40],  This  set 
of  points  is  shown  In  Figure  0.3. 0.0.  Part  of  the  unit  square  Is  outside  the 
convex  hull.  The  points  are  listed  in  Table  3. 

Two  sets  of  data  were  obtained  from  the  literature,  and  one  of  these  was 
scaled  in  one  variable  to  obtain  another.  A  fourth  set  was  used  to  generate 
a  "Cardinal  Function".  The  data  given  in  Table  3,  and  shown  In  Figure  0.3. 0.0. 
was  given  the  following  function  values:  f(x(<,  yk)  ■  0  except 
f(  .1875,  .2625)  =  .2.  Here  .2  was  used  for  visual  purposes  rather  than  1  as 
would  ordinarily  be  done  for  a  true  cardinal  function.  This  gives  some  Infor¬ 
mation  about  the  Influence  of  one  point  on  the  surface  for  moderate  sized 
point  sets.  Of  the  two  sets  of  points  from  the  literature,  one  Is  from  Aklrna 
[1]  and  was  obtained  during  a  study  of  waveform  distortion.  It  Is  repeated 
here  In  Table  5,  and  shown  In  Figure  0.5, 0.0.  The  second  was  obtained  from 
Ferguson  [14]  and  Is  repeated  here  In  Table  6,  and  shown  In  Figure  0.6. 0.0. 

The  same  set  of  data,  but  with  the  y  coordinate  multiplied  by  three  was 
also  used  to  show  effects  of  scaling  only  one  variable,  and  is  shown  In 
Figure  0.7. 0.0.  For  visual  purposes,  the  function  values  given  In  Table  2 
are  actually  .5  more  than  given  by  Ferguson.  As  can  be  seen  from  Figure 
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0,6. 0.0.  the  convex  hull  of  the  data  Is  trapezoidal  shaped.  Since  the 
plotting  routines  expect  function  values  on  a  rectangular  grid  it  was 
decided  to  evaluate  the  Interpolating  surfaces  on  a  rectangle  which 
contained  most  of  the  convex  hull,  but  also  Included  a  mild  amount  of  extrap¬ 
olation.  The  rectangle  was  [2,  18]  x  [-3,  3.4]  on  the  original  data  and 
[2,  18]  x  [-9,  10.2]  on  the  modified  data.  The  convex  hull  of  Akima's 
data  Is  rectangular  and  this  rectangle  was  used  for  evaluating  the  surface 
points. 

The  problem  of  extrapolation  outside  the  convex  hull  has  been  addressed 
by  taking  the  attitude  that  while  It  Is  undesirable  to  have  to  do  so,  It  Is 
likely  possible  to  do  It  In  a  "reasonable"  fashion.  Certainly  In  many 
Instances  (our  cases  mostly  among  them)  one  may  have  better  information  for 
mild  extrapolation  than  for  some  points  within  the  convex  hull.  The  final 
result  Is  that  some  programs  were  modified  to  extrapolate  in  a  "reasonable" 
manner,  some  were  Implemented  that  way  to  begin  with,  and  with  others  the 
problem  does  not  arise.  Basically  only  triangle  based  programs  need  to 
address  the  prob’em,  and  among  those,  only  Lawson's  [33]  program  does  no 
extrapolation.  Points  outside  the  convex  hull  were  omitted  from  the 
deviation  statistics  in  Lawson's  method.  For  the  100  point  data  set,  only 
13  points  of  the  1089  evaluation  points  were  outside  the  convex  hull,  and 
for  the  25  point  data  set  54  points  of  the  1089  evaluation  points  were  outside 
the  convex  hul 1 . 
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1.3.0.  Plot  and  Table  Identification  Scheme 

The  output  of  this  study  consists  In  part  of  a  large  number  of  per¬ 
spective  plots  of  surfaces  and  extensive  tables.  For  ease  In  referencing 
them  they  have  been  gathered  at  the  end  of  the  report,  and  the  entire 
report  has  been  published  In  loose  leaf  form  to  facilitate  reader  compar¬ 
isons  of  corresponding  plots. 

The  plots  have  been  arranged  according  to  a  scheme  Involving  4  numbers, 
each  of  which  identifies  a  particular  aspect  of  the  plot.  The  plot  Ident¬ 
ification  Is  of  the  form  Figure  N^.Ng.N^.N^,  where  the  are  used  to 
Identify  the  characteristics  listed  below. 

N1  -  Type  of  plot 

0  -  plot  of  (x,  y)  point  set 

1  -  Indicates  plot  has  four  3"  plots  as  arranged  In  Figure  1 

2  -  Indicates  plot  has  ..four  3"  plots  as  arranged  in  Figure  2 

3  -  Indicates  plot  has  four  3"  plots  as  arranged  In  Figure  3 

4  -  one  6"  plot  per  page 

N2  -  Indicates  (x,  y)  or  (x,  y,  z)  point  set  used 
0  -  does  not  apply 

1  -  100  points  as  described  In  Section  1.2.2 

2- 33  points  as  described  In  Section  1.2.2 

3- 25  points  as  described  In  Section  1.2.2 

4  -  all  of  1,  2,  3  were  used  as  Indicated  In  Figure  1 

5- 50  points  from  Aklma  [1],  given  In  Table  7 

6- 25  points  from  Ferguson  [13],  given  In  Table  5 

7- 25  points,  obtained  as  Ferguson's  points  with  y  coordinate  x  3. 


< 


Nj  -  Test  surface 

0  -  does  not  apply 

1  -  f^  as  defined  In  Section  1.2.2.,  1  s  1  s  6. 

-  Program  Number,  as  Identified  In  Section  3,  and  given  In  Table  S. 


a 
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Surface 

Hb - 

100  point 
Interpo¬ 
lant 

a 

Cardinal 

Function 

b 

Aklma 

Surface 

c 

33  point 
Interpo¬ 
lant 

d 

25  point 
Interpo¬ 
lant 

c 

Ferguson 

Surface 

d 

Modified 

Ferguson 

Surface 

Figure  1 


Figure  2 


Test 

Surface 

Parameter 
<  Nominal 
Value 

c 

d 

Parameter 

Parameter 

■  Nominal 

>  Nominal 

Value 

Value 

Figure  3  ' 

The  plots  all  Involved  evaluation  of  the  Interpolant  on  a  33  x  33  grid 
of  equally  spaced  points.  Generally  this  grid  Is  fine  enough  so  that  piece- 
wise  linear  plotting  of  the  cross  sections,  which  Is  the  process  used  by 
the  plot  program,  yields  sufficiently  smooth  looking  results.  In  some  Instances 
this  Is  not  really  fine  enough  to  show  the  true  character  of  the  surface,  but 
In  these  f.ases  the  surface  Is  not  a  good  approximation  to  the  test  surface 
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and  the  plot  Is  considered  sufficiently  accurate  to  evaluate  the  visual 
aspect  of  the  surface  anyway. 

Tables  of  comparative  results  are  arranged  l  d  labeled  according  to 
Information  contained  and  test  function  and  data  set  from  which  It  arose, 

If  pertinent.  There  are  several  kinds  of  tables:  (1)  Deviation  tables, 
giving  the  maximum,  mean,  and  root-mean- square  deviations  over  the  set  of 
evaluation  points  used  for  plotting.  These  are  labeled  Table  D.M^Mg,  , 
where  D  Indicates  "deviation  table",  M^  ■  1  means  (x,  y)  data  set  1,  as 
described  for  N2  ■  1,  above;  and  M2  *»  1  Indicates  the  test  function  f^ 
as  described  for  N3  above.  (2)  Timing  tables,  giving  the  execution  times  In 
seconds  on  the  IBM  360/67.  These  times  are  divided  Into  preprocessing  (for 
methods  for  which  there  Is  preprocessing),  evaluation,  and  total.  All  pro¬ 
grams  were  compiled  using  the  Fortran  H  (optimizing)  compiler.  Since  the 
configuration  of  the  machine  Involves  multiprogramming,  these  times  are 
dependent  on  external  factors,  and  may  vary  10$  or  more,  In  either  direction, 
on  otherwise  Identical  runs.  Therefore,  times  are  given  to  two  digits,  the 
second  probably  not  being  significant.  The  tables  are  labeled  Table  T.M, 
where  T  Indicates  "timing  table"  and  M  ■  1  means  for  the  (x,  y)  data  set  as 
described  for  N?  s  1,  above.  (3)  Parameter  variation  tables  give  the  devi¬ 
ations  for  the  nominal  value  of  the  parameter  (for  methods  Involving  a  param¬ 
eter),  and  for  values  larger  and  smaller  than  the  nominal  value.  The  tables 
are  labeled  Table  P.M,  where  P  Indicates  "parameter  variation"  and  M  ■  1 
means  for  the  (x,  y)  data  set  as  described  for  N2  ■  1 ,  above.  (4)  Summary 
table,  Table  S  summarizes  the  pertinent  information  about  all  tested  methods. 
(5)  Two  tables  compact  the  Information  In  the  deviations  tables,  Indicating 
only  which  method  (by  number)  has  the  smallest  deviations  for  each  test 


surface  and  point  set.  The  two  tables  are  for  local  methods  (Table  E.l) 
and  all  methods  (Table  E.2). 


All  tables  listing  results  for  various  methods  are  grouped  into  two  or 
three  separate  groups.  The  first  group  contains  extensively  tested  local 
methods,  the  second  contains  extensively  tested  global  methods,  and  third 
(when  It  appears  )  contains  all  other  methods. 

Certain  Information  In  the  summary  table.  Table  S  needs  additional 
explanation,  in  particular  those  given  letter  grades.  Sensitivity  to 
parameters  Is  a  purely  subjective  score,  based  on  informal  testing  of  the 
scheme.  Included  were  whether  some  value  of  the  parameter  worked  well  for 
a  variety  of  surfaces  for  a  given  set  of  (x,  y)  points,  and  whether  the 
Interpolant  was  stable  with  respect  to  changes  in  the  parameter  from  that 
value.  Complexity  simply  reflects  the  Investigator's  perception  as  to  the 
complexity  of  Ideas  Involved  and  the  ease  of  implementation  Into  a  computer 
program.  Accuracy  Is  again  subjective  and  is  based  on  the  relative  amount 
of  deviation  one  might  expect  from  the  true  surface  for  a  given  method.  Of 
course,  perusal  of  the  deviations  tables  will  reveal  some  methods  do  well 
on  some  surfaces  and  not  so  well  (relatively  speaking)  on  others.  Timing 
Is  relatively  well  defined.  The  first  letter  represents  the  sum  of  the 
evaluation  times,  given  In  Tables  T.l,  T.2,  and  T.3.  Ranges  for  A,  B,  C, 

D,  and  F,  respectively,  are  (0,  7],  (7,  21],  (21,  30],  (30,  50],  and 
(50,  “).  The  second  letter  represents  the  total  time  for  100  data  points 
and  1089  evaluation  points,  the  time  given  in  Table  T.l.  Ranges  are  (0,  4], 
(4,  12],  (12,  20],  (20,  30],  and  (30,  -) 
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2.0.0.  Descriptions  of  Tested  Methods 


For  description  purposes  the  methods  are  classed  into  six  groups: 

(1)  Inverse  distance  weighted  methods,  (2)  Franke's  method,  (3)  Triangle 
based  blending  methods,  (4)  Finite  element  based  methods,  (5)  Foley's 
methods  and  (6)  Nodal  basis  function  methods.  While  there  is  necessarily 
a  blurring  of  distinctions  across  these  group  lines,  they  constitute  fairly 
distinct  Ideas  and  it  Is  convenient  to  group  them  this  way.  In  the  Section 
headings,  the  number  appearing  is  the  number  assigned  to  the  program  Imple¬ 
menting  that  scheme.  This  number  has  no  significance  except  that  It  gives  the 
approximate  order  In  which  the  programs  were  Implemented  or  obtained.  Not 
all  numbers  appear  because  certain  Ideas  were  discarded  as  not  within  the 
context  of  the  study  (In  one  case),  or  as  extremely  deficient  (one  case). 

The  programs  Included  In  the  test  and  a  few  words  describing  it  (also  used 
in  Section  headings  In  this  chapter)  are  given  In  Table  S  to  have  them 
available  for  easy  reference. 
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2.1.0.  Inverse  Distance  Weighted  Methods 


The  original  Inverse  distance  weighted  Interpolation  method  Is  due  to 
Shepard  [SO].  All  methods  of  this  type  which  we  consider  may  be  viewed  as 
generalizations  of  Shepard's  method,  or  variations  of  such  generalizations. 

The  basic  Shepard's  method  Is 

N  N 

(1)  F(x,  y)  =  l  wk(x,  y)fk/  l  wk(x,  y), 

k*l  *  K  k-1  K 

where  v’k(x,  y)  ■  dk’w,  and  typically  y  ■  2,  although  other  values  may  be 
used,  u  may  be  replaced  by  yk  and  could  possibly  bo  different  for  each  k. 
Several  authors  have  considered  various  aspects  of  Shepard's  method  [4],  [5], 
[20],  [49]. 

Shepard's  method  Is  a  global  method,  and  the  original  paper  suggested 
a  scheme  for  localizing  It  by  piecing  together  a  parabolic  segment  with 
dk  In  such  a  way  as  to  obtain  a  wk  which  Is  zero  outside  some  disk,  say  of 
given  radius  R,  centered  at  (xk,  yk),  and  which  Is  still  C1 .  A  simpler  and 
more  natural  scheme  suggested  by  Franke  and  Little  [4,  p.  112]  Is  used  In 
much  of  this  work,  that  Is, 

r<*  -  o+  I2 

(?)  »k<x.  =  1-13"- 

L  K  J 

Shepard's  method  has  an  undesirable  property  for  general  use  In  that  a 
flat  spot  occurs  at  each  data  point.  Use  of  information  about  derivatives, 
either  given  or  generated  from  the  data  was  suggested  by  Shepard,  and  resulted 
in  an  approximation  of  the  form 


(3)  F(x,  y)  *  l  wk(x,  y)[f.  +  (|£)  (x  -  xj  +  (f£)  (y  -  yj]/  l  wk(x,  y). 
k=l  K  K  3X  k  K  3y  k  K  k-1  K 


More  generally,  one  may  consider  approximations  of  the  form 


(4) 


© 


N  N 

F(x,  y)  =  l  w.(x,  y)Ltf(x,  y)/  2  wk(x,  y), 
k*l  K  K  k«l  * 


where  Lkf  Is  an  approximation  to  f  such  that  Lkf(xk,  yk)  ■  This  Is  the 
basis  for  several  of  our  methods.  In  this  context  we  refer  to  the  Lkf  as 
nodal  functions. 


Another  way  In  which  Shepard's  method  can  be  generalized  Is  to  view  the 


method  as  an  Inverse  distance  weighted  least  squares  approximation  to  f(x,  y) 


by  a  constant.  One  can  then  generalize  to  an  approximation  taking  the  form 

(5)  F(x.  y)  ■  P(a0,  a1 . afl;  x.  y), 

where  ag,  . ...  an  are  parameters  chosen  by  taking  them  to  minimize  (for  a 
given  (x,  y ) )  the  expression 


N 

I  [ft.  "  F(aQ,  a-. ...  i 
k«l  *  01 


*n;  xk,  yk)]zwk(x,  y) . 


This  approach  was  taken  by  McLain  [35]  In  evaluating  a  number  of  methods  where 
?  was  taken  as  a  linear  combination  of  low  order  monomials  and  wk(x,  y)  as 
dk  or  exp(-adk  )d,(“  .  McLain  also  considered  some  approximations  where  f 
entered  nonll nearly.  We  have  considered  one  of  McLain's  methods  and  a 


variation  of  another.  All  of  the  methods  of  this  class  may  be  derived  as 
variations  of  the  above  formula  for  P[1 8] . 


We  consider  Shepard's  method  mainly  to  show  how  the  original  method 
performs  In  comparison  with  variations.  The  formula  Is  described  by  Equation 
(1),  but  was  achieved  computationally  as  a  special  case  of  the  Modified 
Shepard's  Method  by  taking  R  very  large. 


2.1.2.  Modified  Shepard's  Method  (7) 

This  variation  is  obtained  by  using  the  weight- function 


wk(x,  y) 


-2 

in  place  of  dk  .  In  general  R  can  be  different  for  each  k,  but  we  have  not 
done  this.  In  order  to  simplify  the  choice  of  R  and  to  remove  effects  of 
scaling  from  the  procedure,  R  Is  actually  computed  from  the  expression 

(6)  R“?,^rD* 


where  D  Is  the  diameter  of  the  point  set  { (xk ,  yk)}  and  Nw  Is  a  new  parameter 
to  be  specified.  Geometrically  Nw  represents  the  anticipated  number  of  points 
which  will  be  in  a  disk  of  radius  R.  Computational  experiments  have  led  to 
a  nominal  (or  default)  value  of  N,.  ■  12.  For  point  sets  of  widely  varying 
density  this  Is  probably  not  an  appropriate  value,  since  the  use  of  constant 
R  for  all  k  assumes  a  somewha,t  uniform  distribution. 
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2.1.3.  Modified  Linear  Shepard'S  Method  (3) 

This  variation  Is  obtained  by  taking  Lfcf  In  Eq.  (4)  to  be  the  Inverse 
distance  weighted  least  squares  approximation  to  the  (xj.  y j »  f4)  by  a  plane 


with  weight  given  by 


and  the  weight  function  Is  that  given  by 


(2).  The  comments  regarding  the  choice  of  R  In  the  previous  method  apply 
here  as  well,  including  the  nominal  (or  default)  value  of  1^  a  12.  The 
coefficients  of  the  plane  are  obtained  In  a  preprocessing  phase. 
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2.1.4.  Modified  Shepard's  Method  Boolean  Sum  Plane  (2) 

Barnh111nand  .Gregory  [6]  have  shown  that  the  operator  PtQ«P  +  Q-PQ 
has  the  interpolation  properties  of  the  operator  P  and  the  precision  of 
operator  Q.  A  suggested  scheme  for  obtaining  polynomial  precision  for 
Shepard's  method  Is  to  use  an  operator  with  linear  precision  for  Q,  while  P 
Is  taken  as  Shepard's  Method.  We  have  used  the  following  scheme.  Let 
L(x,  y)  ■  a(x,  y)x  +  b(x,  y)y  +  c(x,  y)  represent  the  approximation  to  f(x,  y) 

o 

obtained  by  a  least  squares  approximation  with  weight  (R  -  dk)£  for  the  kth 
point,  and  let  Sf  represent  the  Modified  Shepard's  Method  operator  of  Section 
2:1.2.  where  R  above  Is  the  same  as  In  Section  2.1.2.  Then  the  approximation 
Is  F(x,  y)  “  S  •  Lf(x,  y).  Computationally  this  is  achieved  by  S  •  Lf  • 

$(I  -  L)f  +  Lf,  or 

N  N 

(7)  F(x,  y)  »  l  Wk(x,  y)(fk  -  L(xk.  yk))/  l  Wk(x,  y)  +  L(x,  y). 

k«1  K  K  K  K  K 

Thus  the  values  L(xk,  yk)  are  computed  as  a  preprocessing  step,  and  the  two 
terms  in  Eq.  (7)  are  computed  for  the  given  (x,  y). 
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2.1.5.  Modified  McLain  Method  Mg(8) 

McLain's  method  Mg  [38]  Is  of  the  form  given  In  Eq.  (5)  with 

F(Bq,  ay  a2j  x,  y)  -  aQ(x,  y)  +  a1(x,,y)x  +  a2(x,  y)y  and  Inverse  distance 
-2 

weighting  d^  .  We  have  modified  this  by  taking  weighting  given  by 

/ 

,  where  R  Is  again  computed  from  expression  (6)w1th  a  nominal 
value  of  Nw  ■  12. 


'll 


M— 
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2.1.6.  Quadratic  Shepard's  Method  (17) 


This  method  Is  of  the  form  given  by  equation  (4)  where  wk(x,  y)  ■  dk  . 

and  L^x,  y)  is  the  Inverse  distance  weighted  least  square  quadratic  at  the. 

-2 

point  (x^,  yk)  with  weight  dk  .  The  coefficients  of  the  quadratics  are 
obtained  In  a  preprocessing  phase.  This  method  Is  actually  treated  as  a 
special  case  of  the  next  method  with  R  and  Rq  taken  very  large. 
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are  taken  as  r(0  w  x  -,2 

r(Ro  •  v+i 

m  • 

where  R^  bears  the  same  relationship  to  Nq  as  R  to  Nw  In  Eq.  (6),  l.e., 

1  /^o 

Rq  "  2v  Jr  D>  The  11011,1  na1  value  for  was  determined  by  computational 
experiments  and  Is  Nq  ■  18.  If  fewer  than  6 points  lie  in  a  disk  of  radius 
Rq  at  some  (x^,  yk),  the  approximation  Is  taken  to  be  linear.  In  any  case 
nonuniqueness  of  the  nodal  functions  Is  avoided  by  using  the  pseudo-inverse. 


obtaining  the  least  squares  approximation  which  has  minimum  *2  norm  of  the 
coefficients.  The  weight  function  Is  given  by  equation  (2),  and  R  obtained 


from  fL  with  the  nominal  value  of  II.  “  9.  Complete  details  are  given  In 
w  w 


2.1.6.  McLain's  Method  M1(J  (5) 

where 
the  weight 
Here,  In 

D‘ 

where  D  Is  the  diameter  of  the  point  set  { (xk .  yk)>.  This  choice  yields 
a  ■  1  In  McLain's  original  numerical  experiments,  where  McLain  suggested 
a  should  be  something  like  the  usual  distance  to  the  nearest  neighbour. 
Experiments  have  confirmed  that  the  above  a  Is  a  reasonable  choice  In  a 
variety  of  Instances.  This  Is  a  global  method. 


McLain's  Method  M^  [39]  Is  of  the  form  given  by  Eq.  (5) 

F(a0,_  ....  agi  x,  y)  aQ  +  a^  +  a^  +  a3x2  +  a^xy  +  agy2,  and 

2  -2 

function  for  the  approximation  is  taken  to  be  exp(-adk  )dk  . 
order  to  remove  the  effects  of  scaling  we  have  taken 

»  -  L§2 JL, 


2.2.0.  Franke's  Method 


I 


\ 


} 


I 

\ 


The  class  of  methods  [16]  was  Inspired  by  a  short  paper  by  Maude  [37] 
which  generalized  the  Idea  of  deficient  qulntlc  splines  to  several  variables. 
Unfortunately  the  original  Interpolation  function  exhibits  rather  poor 
behavior  and  has  not  even  been  Included  In  our  tests.  The  original  Idea  was 
to  represent  the  Interpolation  function  as 

N  N 

(8)  F(x»  y)  ■  l  Wk(x,  y)Qk(x,  y)/  £  Wk(x,  y), 

k“l  K  K  k-1  K 

where  Qk(x,  y)  Is  the  quadratic  polynomial  Interpolating  f(x,  y)  at  (xk>  yk) 
and  the  five  nearest  neighbors  to  (xk,  y^)  from  the  set  { (xj »  y j ) > .  and 

0  dk  >  Rk  * 

where  Rk  Is  the  distance  between  (xk,  yk)  anu  its  5th  closest  neighbor.  This 

Idea  was  generalized  to  Include  any  Wk(x,  y)  which  have  finite  support  (to 
make  the  method  local)  so  long  as  the  Qk(x,  y)  Interpolate  f(x,  y)  at  all 
(Xj,  7j)  where  Wk(x^,  y^)  *  0.  Use  of  approximations  Qk(x,  y)  In  Hilbert 
spaces,  particularly  In  Sard  spaces,  was  suggested  and  Implemented  [17]*  One 
of  the  chief  advantages  of  this  approach  Is  that  Instead  of  taking  Wk  with 
disks  centered  at  the  (xk,  yk)  as  support  regions,  It  Is  easy  to  use  a  smaller 
number  of  overlapping  rectangles  In  such  a  fashion  that  at  most  four  terms  In 
the  sum  are  nonzero,  and  l  Wk(x,  y)  =  1.  Use  of  rectangles  also  simplifies 
the  problem  of  determining  which  terms  are  nonzero  and  thus  results  In  a 
faster  algorithm.  In  general,  schemes  of  this  sort  are  given  by 
F(x,  y)  -  l  W4(x,  y)QJt(x,  y),  with  l  W£  5  1  and  certain  Interpolation  conditions 
Imposed  on  the  Q£. 


Wk(x,  y)  ■ 
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The  details  of  the  rectangle  selection  process  follow.  As  an  option 
the  user  may  specify  rectangle  boundaries,  however  an  automatic  selection 
process  Is  available  and  Is  assumed  to  be  the  usual  option.  A  parameter 
NPPR  (for  number  of  points  per  rectangle)  Is  specified,  the  suggested  value 
being  NPPR  »  6.  In  the  automatic  case  we  take  nx  *  ny  ■  [2  /N/nPpr  -  j]  and 
then  grid  lines  In  the  x  direction  at  xQ,  x^  ....  xn  +1  are  chosen  so  that 
each  strip  (x^,  x1)x(-»,  »)  contains  approximately  N/(nx  +  1)  points.  A 
similar  partition  Pq,  y^,  . ..,  yn  Is  found  In  the  y  direction.  Now, 
weight  functions  W^x,  y)  with  support  [x^,  x1+1  ]x[pj_., ,  yj+1]  ■  R^ 
are  used  together  with  Q^x,  y)  which  satisfy  Q1J(x|<,  yk)  -  fk' whenever 
(xk,  yk)  e  R^j  to  form  the  Interpolation  function 

(9)  F(x,  y)  »  l  W^{x,  y)Qi(j(x,  y). 

1 1  j 

Here  we  choose  the  so  that  l  si. 

Recently,  some  work  due  to  Junklns,  Jancaltus,  and  coworkers  [31],  [33] 
has  come  to  the  Investigator's  attention.  This  work  involves  the  Idea  of 
weighted  local  approximations  In  a  similar  fashion,  and  was  applied  to  the 
problem  of  terrain  modeling.  In  their  case  the  local  Interpolation  functions 
were  replaced  by  least  squares  approximations  by  polynomials  and  thus  Inter¬ 
polation  was  not  achieved. 


i 
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2.2.1.  Franke's  Method  (Mode  One)  (6) 

7n  this  method  the  Q^j  of  Equation  (9)  are  taken  as  the  optimal  Inter¬ 
polation  function  in  the  Sard  corner  space  .  These  functions  are 
bicubic  spline  functions  and  have  continuous  second  derivatives  except  along 
two  lines  x  ■  a  and  y  ■  b.  By  taking  (a,  b)  outside  of  the  rectangle  R^ 
the  function  Is  then  C  on  R^j.  To  preserve  the  approximation  under 
scaling  (not  necessarily  the  same  In  each  variable)  the  optimal  approximation 
Is  computed  after  R,y  Is  transformed  to  [0,  1]  .  At  least  three  Interpolation 
points  are  used,  nearest  points  (In  the  norm)  being  added  If  necessary. 

To  preserve  the  continuity  of  the  second  derivative  It  Is  necessary  to  take 
W^j  with  continuous  second  derivatives.  Thus  the  choice  of 

Wij(x,  y)  *  Vi(x)Uj(y),  where 
1  ,  x  <  x. 


V,(x) 


X  -  X, 


H  (-  ,  x,  S  X  <  Xr 


x2  -  X, 


X  *  x. 


v1(x) 


v„  <x) 

X 


,  X  <  X 


1-1 


1  *  V,j_.|  (x) ,  x1_1  S  x  <  x1 


h5<“ 
0 


X  -  X, 


-).  Xi  S  X  <  X, 


xi+l  ‘  X1 


1H 


X  2  X 


1+1 


for  1  *  2 ,  . . . ,  n  -  1 ,  and 

0  •  x  <  Vi 

1  -  Vn  -l(x),  X  -1  <  x  x 
'x  x^  x 

1  ,  X  >  X  , 

X 
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is  made.  where  Hg(s)  ■  l  -  s3(6s2  -  15s  +  10)  is  the  HSrmlte  quint 1c 
satisfying  H,(0)  -  1,  H£<0)  -  Hg(0>  -  Hs(l)  .  HJ(1)  »  H|(l)  -  0. 


2.2.2.  Pranke's  Method  (Mode  Three)  (1) 

Because  the  optimal  approximations  In  Bp2  gj  have  no  polynomial 
precision,  another  choice  for  local  approximating  functions  is  available. 
In  this  case  the  j  are  taken  to  be  the  optimal  approximation  In 
boolean  sum  the  least  squares  (unweighted)  plane  fit  to  all  data  points  In 
R^j.  Since  the  latter  process  has  linear  precision,  so  does  the  overall 
approximation.  The  process  Is  Implemented  as 

B#Lf  *  B(I  -  L)f  +  Lf,  where 

B  Is  the  optimal  approximation  and  L  Is  the  least  squares  plane  fit. 

The  choice  of  rectangles  and  weight  functions  Is  Identical  to  that  of 
the  previous  section. 
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The  elegant  theory  and  excellent  fitting  characteristics  of  the  thin 
plate  approximations  given  by  Duchon  [1*1],  (see  Schumaker  [47]  and  Section 
2.5.4.)  lead  to  their  consideration  as  local  approximations  In  the  basic 
method  given  by  Equation  (9).  Several  other  modifications  were  incorporated 
as  well. 

The  suggested  number  of  subrectangles  remains  the  same.  However  the 
selection  of  grid  lines  x^  and  yj  Is  done  in  a  way  which  preserves  symmetry 
under  reflections  and  also  will  result  in  a  symmetric  Interpolant  If  the 
data  Itself  Is  symmetric. 

In  selecting  points  to  be  Interpolated  by  the  a  slightly  larger 
rectangle  than  Is  considered  by  Including  all  points  In  the  rectangle 
[-.1125,  1.1125]2  after  has  been  transformed  to  [0,  l]2,  This 
rectangle  has  area  approximately  50%  larger  than  unity  and  Interpolation  on 
the  larger  set  of  points  tends  to  make  the  transition  between  regions  somewhat 
smoother.  This  choice  was  made  on  the  basis  of  computational  experience.  Agal 
at  least  three  points  must  be  Interpolated  and  the  nearest  points  (In  the 
norm)  are  added  If  necessary. 

Experience  has  shown  that  many  C1  surfaces  appear  to  be  smoother  than  C1 
In  that  second  derivative  jumps  are  apparently  small.  While  the  thin  plate 
approximations  have  discontinuous  second  derivatives  at  the  data  points,  the 
former  reason  Is  the  primary  one  for  using  H3(s)  =  1  -  s  (3  -  2s)  In  place  of 
H5(s)  In  the  definition  of  the  W^  l-’or  this  method. 

The  local  approximations  have  the  form 


where  Is  the  set  of  Indices  k  for  which  (xk,  yk>  fk)  is  a  point  to 
be  Interpolated  by  Q^j.  See  Section  2,5.4.  for  a  further  discussion  of 
thin  plate  splines. 


2,3.0.  Triangle  Based  B'lendl rtg  Methods 

These  methods  are  conceptually  the  same  as  given  by  Equation  (4), 
but  a  significant  difference  Is  that  the  weight  functions  are  based  on  a 
trlangulatlon  of  the  convex  hull  of  the  point  set  {(xk,  y^}}.  Several  such 
schemes  have  been  proposed,  e.g.  [9],  [18],  [19],  and  [40],  One  of  those 
considered  here  Is  the  one  described  In  [18].  , 

Assume  a  trlangulatlon  of  the  convex  hull,  and  suppose  (x,  y)  e  T.j where 
T1jk  Is  the  triangle  with  vertices  (x^,  y^),  (Xj ,  y j ) ,  and  (xk,  yk).  We  then 


(10)  F(x,  y)  »  W^x,  yjQ^x,  y)  +  W^(x,  y)Qj(x,  y)  +  Wk(x,  y)Qk(x,  y) 
where  the  weight  functions  are  finite  element  "shape"  functions  satisfying 

VV  “  5mt  end  V  V  “  f4  for  m,  i  ■  1,  j,  k.  In  all  previously 
referenced  methods  the  weight  functions  may  be  vlewod  as  nine  parameter  cubic 
shape  functions  with  a  rational  correction  to  obtain  normal  derivatives  equal 
to  zero,  and  hence  a  C1  approximation  overall.  There  are  many  ways  to  obtain 

*  i 

such  correction  terms,  all  of  which  appear  to  lead  to  the  posslbllty  of  negative 
values  being  taken  on  by  one  of  the  weight  functions  If  the  triangle  is  very 
obtuse.  This  Is  probably  not  serious,  although  one  has  no  control  over  the 
shape  of  the  triangles  In  the  sense  that  very  obtuse  angles  cannot  be  avoided 
In  some  instances.  The  weight  functions  used  here  are  obtained  from  a  minimum 
norm  problem  [43].  Let  ,  bj,  hk  be  the  barycentrlc  coordinates  of  (x,  y) 

In  V  and  let  Hy  z j ,  and  &k  be  the  lengths  of  the  sides  opposite  vertices 
1,  j,  and  k,  respectively.  Then  the  weight  function  is  given  by 


Wk(x,  y)  «  bk(3  -  2bk)  +  Sb^b^a^  +  ak1] 


W  +  M 
n  -  bv  m  - 


2  2 
V  £1 


-39- 


and  the  others  being  obtained  by  a  cyclic  permutation  of  the  Indices. 

While  the  basic  method  Is  defined  only  on  the  convex  hull  of  the  point 
set,  It  Is  easily  extended  to  be  a  globally  defined  function  by  the  following 
Idea,  The  exterior  of  the  convex  hull  Is  divided  Into  semi Infinite  rectangles, 
shown  In  Figure  4,  by  iconstructlng  perpendiculars  to  the  exterior  edges  of 
the  convex  hull  at  each  exterior  vertex. 


Figure  4 


To  extend  the  definition  of  Equation  (10)  outside  the  convex  hull  we 
proceed  as  follows.  For  a  point  In  an  exterior  triangle,  such  as  (x,  y), 


we  take  F(x,  y)  ■  Qj(x,  y).  For  a  point  In  an  exterior  rectangle,  such  as 
•  »  *  » 

(x,  y),  let  p  be  the  projection  of  (x,  y)  onto  side  1j,  and  let  (bj,  bj,  0) 


be  the  barycentrlc  coordinates  of  p  In  T^. 


Then  we  take 


F(x,  y)  ■  h3(b1)Q1(x,  y)  +  h3(bj)Qj(x,  y), 


p 

where  h3(s)  ■  s  (3  -  2s).  These  extensions  yield  a  globally  defined  approx- 
Imatlon  which  Is  C^. 
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2.3.1.  Nlelson-Franke  Linear  Triangle  Method  (12) 

This  method  uses  the  Inverse  distance  weighted  least  squares  plane 
Qk(x,y)  which  Is  also  used  In  the  Modified  Linear  Shepard's  Method.  See 
Section  2.1.3.  for  details. 

Another  Idea  was  Investigated  for  determining  the  slopes  to  be  used  In 
the  planar  fit,  but  was  abandoned  as  being  a  poor  Idea.  Since  the  Idea  has 
been  mentioned  by  a  number  of  persons,  It  Is  discussed  here.  At  a  vertex  k, 
determine  the  slopes  a^  and  b^  of  the  plane  a^x  +  b^y  +  c^j  through  the 
points  (x.|,  y^ ,  f^),  (Xj»  yj.  fj)»  ®^d  (*k*  yk.  fk)»  where  ^jk  Is  a  triangle 
In  the  trl angulation.  The  nodal  function  Is  taken  as  Qk(x,  y)  ■  Ak(x  -  xk) 

+  Bk(y  -  yk)  +  fk,  where  Ak  and  Bk  are  the  average  values  of  a^  and  b^, 
respectively.  We  can  think  of  this  as  taking  the  nodal  function  at  each  vertex 
as  the  average  plane  from  the  piecewise  linear  approximation  on  the  trl angulation. 
This  falls  because  of  the  possible  existence  of  long  thin  triangles.  This  Is 
especially  crucial  when  the  triangle  Is  very  obtuse,  and  the  plane  through  the 
three  points  may  have  very  large  gradients  because  the  three  vertices  lie 
nearly  on  a  straight  line  while  the  three  points  on  the  surface  do  not. 


2.3.2.  Nlelson-Franke  Quadratic  Triangle  Method  (13) 

This  method  uses  the  Inverse  distance  weighted  quadratic  Qk(x,  y)  which 
Is  also  used  In  the  Modified  Quadratic  Shepard's  Method.  See  Section  2.1.7. 
for  details.  A  complete  description  Is  given  In  [18]. 
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2.4.0.  Finite  Element  Based  Methods 


These  methods  are  based  on  the  concept  of  using  C1  finite  element 
functions  on  a  trlangulatlon  of  the  convex  hull  of  the  point  set.  This 
requires  a  scheme  for  estimating  some  derivatives  (which  ones  depend  on 
the  element  used  by  the  method)  at  the  data  points.  Our  test  results  Indicate 
that  acpuratq  estimates  of  the  derivatives  are  very  Important  and  have  a 
pronounced  effect  on  the  visual  aspects  of  the  surface,  particularly,  but  r 
also  the  accuracy. 


2.4.1.  Aklma 1 s  Method  (6) 

Aklma 's  method  [1]  uses  the  18  parameter  qulntlc  element  and  thus 
requires  both  first  and  second  partial  derivatives  at  each  data  point. 

The  scheme  used  Is  as  follows.  The  user  specifies  a  parameter,  nd.  Let 
Pk  “  Uk*  yk.;>zk)  ,*nd  form  the  sum  V  «  £  ±  x  ,  where  (1,  j),  1  +  j 
ranges  over  the  nd  closest  points  (x^,  y^)  and  (Xj,  yj)  to  (x^ •  yk),  and 
where  the  sign  Is  chosen  so  that  the  z  component  of  each  cross  product  is 
positive.  The  first  partial  derivatives  are  taken  to  be  those  of  the  plane 
normal  to  V.  The  second  derivatives  are  obtained  by  applying  the  same 
process  to  the  derived  data.  The  cross  partial  Is  taken  as  the  average  of 
the  two  so  obtained.  Aklma  suggests  nd  *  3  or  4  as  appropriate,  but  we  have 
found  nd  ■  6  generally  works  better. 

Extrapolation  outside  the  convex  hull  Is  achieved  by  construction  of  an 
appropriate  polynomial  In  the  exterior  rectangle  or  triangular  regions  given 
In  Figure  4,  and  C1  continuity  Is  maintained.  In  a  triangular  region,  the 
conditions  at  the  vertex  determine  a  unique  bivariate  quadratic.  For  the 
rectangular  region  the  conditions  at  the  two  vertices  determine  a  unique 
polynomial  of  degree  two  In  the  direction  normal  to  the  boundary  segment, 
matching  the  quadratic  in  the  adjacent  triangular  region,  and  of  degree  5 
in  the  tangential  direction,  matching  the  value  and  first  two  derivatives 
across  the  boundary  from  rectangular  to  triangular  region. 
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2.4.3.  Aklma's  Method  -  Modification  Two  (11) 

It  was  felt  that  perhaps  an  Inverse  distance  weighting  of  the  cross 
products  might  be  desirable,  so  this  scheme  formed  the  vector  sum 

'  j.E-ifZL 


il 


2 

All  other  aspects  of  the  program  were  maintained, 


2.4.2.  Akima ' s  Method  -  Modification  One  (10) 

It  Is  easily  observed  that  while  Aklma's  scheme  for  estimating  derlv 
atlves  puts  less  weight  on  nearly  colllnear  points,  which  seems  desirable 
It  also  puts  more  weight  on  distant  points,  which  does  not  seem  to  be 
desirable.  To  renove  the  distance  weighting,  the  scheme  was  modified  to 
form  cross  products  of  unit  vectors  In  the  same  directions  as  before,  1,e 
the  sum  >1 


was  formed.  All  other  aspects  of  the  program  were  maintained. 


/ 


2.4.3.  Aklma *s  Method  -  Modification  Two  HP 

It  was  felt  that  perhaps  an  Inverse  distance  weighting  of  the  cross 
products  might  be  deslrablei  so  this  scheme  formed  the  vector  sum 

2  ■  i\w' 

All  other  aspects  of  the  program  were  maintained. 


,  „  I, i,.l  :!«J>V«l3T 


ij  1,1 


2.4.4.  Aklma's  Method  -  Modification  Three  CIS) 

This  modification  Incorporated  use  of  the  Inverse  distance  weighted 
quadratic  least  squares  polynomial  fit  used  In  the  modified  Quadratic  Shepard 
method  described  In  Section  2.1.7.  The  required  derivatives  were  then  taken 
from  the  quadratic  nodal  function  computed  In  this  manner.  All  other  aspects 
of  the  program  were  maintained. 


2.4.5.  Nielson's  Minimum  Norm  Network  (19) 

This  scheme  [43]  uses  a  cubic  element  with  a  rational  correction  term 

to  obtain  C1  continuity.  One  of  the  basis  functions  Is  that  described  In 

Section  2.3.0.  Only  first  derivatives  are  required.  These  are  obtained  by 

h2  - 

minimizing  the  value  of  /[-**■  F(x(s),  y(s))]  ds,  where  the  Integral  Is  over 

ds6 

the  entire  network  of  edges  In  the  trlangulatlon.  We  note  that  the  Inter¬ 
polation  function  Is  a  univariate  Hermlte  cubic  polynomial  along  each  edge. 
This  Is  a  global  method, 

The  original  scheme  Is  not  able  to  extrapolate  outside  the  convex  hull, 
but  the  following  Idea  was  Incorporated  to  achieve  extrapolation.  In  a 
triangular  exterior  region  as  In  Figure  4,  the  function  Is  taken  to  be  the 
linear  function  determined  by  tihe  value  and  slopes  at  the  vertex.  In  the 
rectangular  region  the  function  Is  extended  by  extrapolating  from  the  pro- 
jectlon  point,  p,  with  t(je  given  slope  and  value  along  a  straight  line.  The 
resulting  surface  Is  only  C°  across  exterior  rectangular  to  triangular  bound¬ 
aries,  but  for  mild  extrapolation  this  will  likely  not  be  noticeable.  An 
appropriate  rational  correction  could  probably  be  made  In  the  triangular  area 
to  achieve  C1  continuity. 
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2.4.6.  Lawson's  Method  (28) 

This  method  is  somewhat  similar  to  Aklma's  in  philosophy  except  for  the 
particular  finite  element  used  and  the  manner  of  estimating  derivatives  [34]. 

A 

The  Clough-Tocher  cubic  element,  which  requires  first  derivatives  at  the 
vertices,  is  used.  The  derivatives  are  estimated  by  fitting  an  inverse  dis¬ 
tance  weighted  quadratic  at  each  vertex.  The  program  Is  presently  not  Set  up 
to  extrapolate  outside  the  convex  hull  of  the  point  set,  although  a  scheme 
for  extrapolation  similar  to  that  used  in  Pile! son's  program  couM  be  incor¬ 
porated.  Time  did  not  permit  this,  however. 


2.5.0.  Foley's  Methods 

Foley's  methods  [15]  Involve  severe!  Ideas.  The  use  of  a  generalized 
Newton  type  Interpolant  Is  Involved  In  them  prominently  and  this  Idea  Is 
discussed  In  Section  2.5.1. 

Another  Idea  which  Is  exploited  successfully  Is  that  of  using  one 
Interpolant  to  generate  a  grid  of  points  on  which  product  type  approximations 
can  be  constructed.  The  product  approximation  will  not,  In  general.  Inter¬ 
polate  the  given  data.  Hence  a  correction  based  on  the  original  approximation 
Is  made  to  the  error.  This  process  Is  termed  a  "delta  sum"  by  Foley,  written 
PAQ,  defined  by  PaQ  *  P8QP,  and  Implemented  as  (PAQ)f  *  P(I  -  QP)f  +  QPf. 

The  Idea  has  greater  generality  than  considered  by  Foley,  but  the  appli¬ 
cation  of  It  seems  to  be  the  appropriate  one.  He  considers  cases  where  the 
product  type  approximation  (taking  the  part  of  Q)  is  either  the  bivariate 
product  Bernstein  .polynomial  or| the  bivariate  product  natural  bicubic  spline. 
The  first  Interpolant  (taking  the  part  of  P)  Is  taken  as  either  the  generalized 
Newton  Interpolant,  or  a  form  of  Shepard's  method.  The  delta  sum  Idea  Is 
applied  In  Iterated  form  for  two  methods. 
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2.5.1.  Generalized  Newton  Interpol  ant  (25) 

The  generalized  Newton  Interpol ant  as  considered  by  Foley  takes  the 
form  F(x,  y)  ■  TN(x„  y).  where  T^x,  y)  ■  f,,  Tk(x,  y)  ■  Tk_.j(x,  y)  + 

f^  "  i(Xl»  yJ.  '  >■■  •  k-1 

Wk(k,  y)  -*"  wk(x~k,  ykT~  *  k  "  2*  ••••  N*  where  'Vx»  "  j1T1  Lj^*  w1th 

d^  V(x  -  xk)*  +  (y  -  yk)Z  .  and  Lj (0)  «  0,  Lj ( t)  *  0  if  t  *  0.  Foley  takes 

•  t2 

Lj(t)  ■  -TT— — -k  where  the  value  of  r.  was  obtained  In  an  ad  hoc  fashion  as 
J  r  +  r/  3 


d 


rj  •  [25dj1  +  13dj2  +  7dj3  +  3d^].  Here  dj^  represents  the  distance  to 

the  1—  nearest  neighbor  to  (Xj,  yj)  In  the  set  of  points  {(xk,  yk)}. 

Unfortunately,  unlike  the  univariate  Newton  polynomial,  this  function 
depends  on  the  ordering  of  the  data  points.  A  number  of  experiments  lead  Foley 
to  two  ordering  schemes.  Let  (x,  y)  be  the  centroid  of  the  set  {(xk,  yk)}»  l.e. 


(x,  y)  “  J-  yxk,  yk).  Let  ak  »  (x  -  xfc)c  +  (y  -  yk)£,  and  arrange  the  points 

In  Increasing  order  of  ak-  This  Is  called  "inulde  out"  ordering,  and  the 
opposite  order  Is  called  "outside  In".  The  two  Interpolants  based  on  these 
orderings  are  called  TI0(f)  and  TOI(f),  respectively,  by  Foley.  Since  each 
appears  to  work  better  In  the  region  from  which  the  final  points  come,  l.e., 

TIO  Is  better  In  the  outer  regions,  while  TOI  is  better  In  the  central  regions, 
a  blending  of  the  two  is  taken  as  the  final  Interpolant.  The  weighting  function 
Is  given  by 

BL(x,  y)  -  l  ■  k.:..«4±i g-ilt! - ,  , 

'  (x  -  xr  +  (y  -  y)  +  d* 

pi  —  2  —  2 

where  D  ■  £  max  [(xk  -  x)  +  (yk  -  y)  ].  The  final  Interpolant,  TF,  Is  then 
given  by 

F(x,  y)  -  BL(x,  y)TIO(f)  +  (1  -  BL(x,  y))T0I(f). 


? 


t 
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2.5.2.  TF  Delta  Sum  Bernstein  Polynomial  (26) 

The  only  additional  Information  In  the  Implementation  of  this  scheme  is 

the  region  on  which  the  data  for  the  Bernstein  polynomial  Is  to  come  (a 

square  In  the  bivariate  case)  and  the  degree  qf  the  polynomial.  Foley  takes 

the  square  to  be  the  smallest  square  containing  the  set  {(xk,  yk)>,  although 

notes  It  might  be  better  to  do  otherwise  In  some  circumstances.  The  degree 

Is  taken  to  be  10,  although  this  means  a  grid  of  11  x  11  points  Is  used  for 

the  Bernstein  approximation.  We  have  followed  Foley,  but  It  might  be  more 

0 

reasonable  to  use  an  m  x  m  grid  where  m  «  N,  as  Is  done  In  the  next  section. 
Let  the  Bernstein  polynomial  for  g(x,  y)  on  (a,  b)  x  (c,  d)  be  denoted  by 
BRN(g).  Then  the  TF  delta  sum  Bernstein  polynomial  approximation  Is  given  by 
TFABRN(f)  ■  TF(f  -  BRN(TF(f)})  +  BRN(TF(f)). 
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The  TF  Interpol ant  Is  also  applied  In  conjunction  with  the  natural 
bicubic  spline.  We  again  need  to  give  a  grid  of  points  on  which  to  compute 
the,  bicubic  spline.  The  selection  routine  supplied  by  Foley  (but  noi  used 
for  the  examples  In  his  thesis)  Is  as  follows.  Let  m  ■  [rlf]  (here  [*]  denotes 
the  Integer  part)  ,  and  k  j].  Sort  the  x-coordlnates  so  that 

X|  s  i<2  <  Xj  i  .  .  .  s  xN>  Then  the  x-grld  lines  are  given  by 

kj 

x.4  I  x<  for  j  ■  1 ,  . . . »  m  -  1 , 

i  k  1«(j-1)k+l  1 

~  1  M 

The  y-grld  lines  are  formed  In  dual  fashion. 

We  now  consider  applying  the  delta  sum  in  Iterated  fashion  to  obtain  a 
sequence  of  operators  GQ,  G^,  ...  .  Let  B  represent  the  natural  bicubic  spline 
operator  on  the  above  grid.  Then  let  GQ  »  TF,  and  successively  form  operators 

Gn+1  “  TF®B6n*  n  "  °'  1 . 

The  calculation  can  then  be  organized  as  follows:  Compute  the  current 
approximation  at  the  grid  points;  construct  the  natural  bicubic  spline  Inter- 
polant  for  the  grid;  correct  the  spline  interpolant  for  the  grid  to  obtain 
Interpolation  at  the  scattered  points  by  adding  In  the  TF  Interpolant  for  the 
'error.  Computationally  this  all  amounts  to  writing  W  ,s  BV  +  TF"  -  BGn>f' 
Under  certain  conditions  the  Iteration  may  converge,  and  can  converge  to  a 
bicubic  spline  function  which  Interpolates  the  original  data.  In  other  Instances 
the  Iteration  appears  to  diverge,  unfortunately.  We  have  taken  3  as  the  nominal 
number  of  Iterations  to  be  used. 


Method  Delta  Sum  Bicubic  Spline  (31) 
The  basic  Idea  of  the  previous  section  was  also  applied  using  a  modified 
Shepard's  method  In  place  of  the  TF  Interpolant.  The  stable.  If  somewhat 
undesirable,  behavior  of  Shepard's  method  would  appear  to  be  well  sultdd  for 
this  use. 

The  basic  modification  to  Shepard's  method  was  one  to  force  a  diminished 

3. 


region  of  Influence  on  the  points,  taking  the  weights  to  be 


2  „  ttD*  n2 

r  "  i r , D 


d/(d/  +  r£) 


,  where 


2  „  2  uk  'uk  T  '  ' 

max  [(xj  -  x)  +  (y,j  -  y)  3.  and  (x,  y)  Is  the  centroid  of  the 


set  { (Xj .  yj)},  as  In  Section  2.5.1,  The  modified  Shepard's  method  used  here 


Is  of  the  form 


F(x>  y)  •  j,  j,ppxbr  "here 


Pk(x,  y) 


dk  (dk  +  r  > 


The  Iterated  delta  sum  Interpolant  Is  then  formed  In  exactly  the  same 
manner  as  In  the  previous  section,  with  the  modified  Shepard  operator,  SM, 

replacing  TF.  Thus  we  have  GQ  «  SM,  Vi  *  SM#BGn,  n  ■  0,  1 .  We  have 

again  taken  3  Iterations  as  the  nominal  value,  and  comments  regarding  conver¬ 
gence/nonconvergence  of  the  Iteration  of  the  previous  section  apply  here. 


2.6.0,  Global  Basis  Function  Type  Methods 

These  methods  can  be  characterized  by  the  following  Idea.  For  each 
(xk,  yk)  simply  choose  some  function  Gk(x,  y),  and  then  determine  coefficients 
Ak  so  that  F(x,  y)  ■  £  AkGk(x,  y)  Interpolates  the  data.  Schemes  which  work 
are  not  so  simple  In  that  appropriate  choices  of  function  Gk  are  not  partic¬ 
ularly  easy  to  make.  Even  If  the  functions  Gk  have  only  local  support  the 
methods  are  global  and  further  they  require  solution  of  a  system  of  N  linear 
equations.  In  all  Instances  we  consider,  the  systems  have  a  symmetric  coefficient 
matrix  (G^ (xj ,  yj))»  but  this  need  not  be  the  case.  Usually  the  Gk  are  really 
functions  of  one  variable,  dk  ■  /(x  -  xk)z  +  (y  -  yk)2.  While  It  seems  that 
functions  Gk  which  diminish  as  one  moves  away  from  the  point  (xk,  yk)  would  be 
best,  this  has  not  been  borne  out  computationally.  Numerous  colleagues  have 
suggested  (among  others)  B  splines,  Gaussian  distributions,  and  other  basis 
functions  which  seem  to  have  an  at  best  shaky  mathematical  justification.  We 
Investigated  several  methods  of  this  type  and  have  found  them  to  work  better 
than  expected.  They  are,  as  mentioned,  global  methods. 


2.6.1.  Rotated  Gaussian  (201 

This  scheme  Is  mentioned  by  Bolandl,  etal.  [7],  Arthur  [3],  and  more 

recently  was  rediscovered  with  a  slight  variation  [48].  It  consists  of  using 

2  2  :l 
■  exp  (-dk  /R)»  where  R  Is  taken  to  be  a  constant  for  all  k.  The  method 

Is  quite  sensitive  to  the  choice  of  R  and  yields  poor  results  with  ease,  but 

will  yield  quite  good  results  for  an  appropriate  value  of  R.  We  used  a 

nominal  value  of  R  ,  where  D  is  the  diameter  of  the  point  set.  The 

2^T 

factor  represents  (approximately)  the  radius  of  a  disk  In  which  one 


could  anticipate  finding  one  point  of  the  set  and  In  some  sense  Is  proportional 
to  the  mean  distance  to  the  nearest  neighbor. 


2.6.2. 


This  method  has  been  used  extensively  by  Hardy  and  his  coworkers  [22-28] t 
In  geographic  and  delated  applications.  The  basis  function  used  Is  the  upper 
hyperboloid  Gk  “  ( (x  -  xk)2  +  (y  -  yk)2  +  r2)1//2,  where  r  Is  a  parameter  which 
determines  the  semi-axis  of  the  hyperbola.  Hardy  [26]  Indicates  that  the  best 
value  for  r  Is  approximately  .815d,  where  d  Is  approximately  the  mean  distance 
to  the  nearest  neighbor.  We  have  not  verified  this  and  have  used  a  nominal 
value  of  r  ■  2.5R,  where  R  Is  the  radius  of  the  disk  which  could  be  anticipated 
to  contain  one  point.  The  actual  parameter  used  by  the  program  Is  NPPR  and 
the  value  of  r  Is  computed  from  r  *  R,  R  ■  £  D/rfl  where  D  Is  the  diameter 
of  the  point  set  { (xk ,  yk)}»  and  a  nominal  value  of  25  Is  assumed  for  NPPR. 

We  observe  better  results  are  generally  obtained  with  larger  r,  but  this  also 
leads  to  poorer  conditioning  of  the  coefficient  matrix  y j ) ) *  and  we 

have  compromised  on  the  above  value  which  corresponds  to  approximately  1.6d. 
Because  of  the  scattered  nature  of  the  data  this  may  vary. 


rnnlaMWlMval  tV  (Hi  wri»if  t  iimiMii 


2.6.3.  Hardv's  Reciprocal  Multi Quadric  (27) 

2 

For  this  method  the  reciprocal  hyperboloid  G^x,  y)  ■  ((x  -  xk)  + 

2  2  1/2 

(y  -  yk)  +  r  )  '  Is  used.  The  value  of  r  used  Is  the  same  as  that  for 
the  previous  method. 
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2,6.4.  Duchon's  Radial  Cubic  Method  (22) 

Duchon  [12]  has  derived  this  method  as  the  optimal  solution  In  a  certain 

Hilbert  space  via  construction  of  the  reproducing  kernel  (see  [49]  for  some 

details).  For  practical  purposes  the  user  must  solve  a  system  of  the  form 
n 

l  Aj(aj(x1y1)  +  axi  +  by1  +  c  -  f  1  *  1  ■  1»  ....  n 
J*1 

l  AjXj  ■  0,  l  A1y1  «  0.  I  A,  »  0* 

J#1  3  J  j-1  J  J  j-1  3 

where  Gk(x,  y)  ■  ((x  -  xk)2  +  (y  -  yk)2)3//2.  Gk  Is  seen  to  be  of  the  form 

3  4 

dk  ,  where  dk  Is  distance  from  (xk,  yk)*  hence  rr\y  name  for  It. 


m 


IVi 
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2.6.5.  Duchon's  Thin  Plate  Function  (23) 

This  method  Is  similar  to  the  previous  one  In  that  It  Is  the  optimal 

solution  In  some  Hilbert  space.  This  one  Is  particularly  Interesting  since 

overall  Interpolating  functions  In  the  Hilbert  space  It  minimizes  the  thin 

2  22  2  2 

plate  functional  //  (M)2  +  2(~—r)  +  (Mr)  .  The  form  of  the  solution 
r2  sx  axay  sy* 

had  been  previously  given  by  Harder  and  Desmarals  [21].  The  method  is  also 
discussed  by  Melnguet  [41]  and  its  fitting  properties  In  connection  with 
smoothing  has  been  Investigated  by  Wahba  [52]. 

The  system  of  equations  Is  Identical  In  structure  to  that  of  the  previous 
method  except  that  Gk(x,  y)  «  dk  log  dk?  where  again  dk  Is  the  distance  from 

<v  *k>- 
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2.6.6.  Rotated  B-Sollnes  (29) 

This  Idea  has  been  suggested  by  various  persons.  We  took  the  B-splIne 
based  on  equally  spaced  points,  with  knots  at  ±R,  ±R/2,  and  0,  where  R  was 
chosen  In  a  manner  similar  to  other  schemes  previously  described.  The 
nominal  value  of  R  was  taken  as  where  D  Is  the  diameter  of  the 


point  set,  and  again  the  factor  —  represents  the  radius  o;f  a  disk  In  which 
one  could  anticipate  finding  one  point  In  the  set.  The  basis  function  was  . 
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3.0.0.  Results 

The  results  of  the  study  are  discussed  In  this  section  In  the  same 

.  ■  ■  ■  j  i  1 

sequence  as  the  methods  are  discussed  In  Section  2,  with  corresponding 
subsection  numbers  relating  to  classes  of  methods.  We  confine  our  comments 
here  to  absolute  merits,  of  each  method,  Insofar  as  possible,  although  It  Is 
almost  certain  we  are  "grading  on  the  curve".  Merits  of  some  methods  (the 
better  ones  or  more  available  ones.  In  our  opinion)  versus  other  methods  are 
discussed  In  Section  4. 

It  Is  hardly  possible  to  discuss  the  performance  of  each  method  on  each 
surface  In  the  detail  which  would  be  desirable  from  a  completeness  point  of 
view,  If  not  the  writer's  (and  perhaps  not  the  reader's  either),  tn  order  to 
point  out  pertinent  behavior  which  one  can  look  for  in  various  methods,  It 
was  decided  to  discuss  In  some  detail  all  plots  of  surfaces  for  one  method. 

In  addition,  there  are  a  number  of  comments  about  the  test  surfaces  relative 
to  the  data  sets  which  apply  to  all  methods.  The  method  dlcussed  Is  neither 
the  best,  nor  the  poorest,  but  simply  a  method  which  Illustrates  some  of  the 
behavior  one  can  watch  for.  The  method  chosen  for  discussion  Is  program  #30, 
Foley's  Iterated  generalized  Newton  delta  sum  bicubic  spline.  We  discuss 
each  figure  separately. 

Figure  1 .4.1.30.  Part  a  Is  the  test  surface.  Part  b  Is  the  Interpol  ant 

2  2 

based  on  100  points.  The  peak  of  the  test  surface  occurs  near  (|,  ^)  and 
Inspection  of  the  set  of  points  In  Figure  0.1. 0.0.  reveals  a  relatively  large 
gap  In  that  vicinity.  Thus  the  peak  value  of  the  surface  Is  not  well  defined. 
Part  b  shows  the  left  rear  of  the  peak  to  be  the  poorest  portion  of  the  peak 
definition.  At  the  left  front  ofthe  surface  an  undesirable  flattening  of  the 
surface  occurs.  Cross  sections  throughout  the  Interior  appear  to  be  quite 
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smooth.  At  the  rear  of  the  dip  there  is  an  extended  depressed  region. 

Reference  to  Figure  0.1. 0.0.  again  reveals  that  there  Is  a  relatively  large 
gap  In  the  points  which  causes  the  extent  of  the  dip  to  be  poorly  defined. 

The  dip  Is  near  (y,  ^).  Figure  4.1.1.30  Is  a  larger  plot  of  the  surface 
In  Figure  1.4.1.30b.  Part  c  shows  several  111  defined  portions  of  the  surface 
corresponding  to  regions  with  gaps  In  the  data  set,  shown  In  Figure  0.2. 0.0. 

In  particular,  near  the  rear  edge  of  the  surface  behind  the  peak,  and  to  the 
right.  The  dip  Is  completely  missed  because  of  a  lack  of  data  to  define  It. 

The  surface  basically  appears  reasonable,  except  possibly  for  the  behavior 
at  the  rear  edge  near  the  right.  A  larger  plot  of  the  data  Is  shown  in  Figure 
4.2.1.30.  The  surface  in  part  d  shows  basically  appropriate  behavior.  The 
peak  Is  reasonably  well  defined,  although  slightly  low.  Again,  no  point  Is 
on  the  top  of  the  peak.  The  dip  Is  somewhat  defined,  but  spread  out  because 
of  a  lack  of  nearby  points  to  pull  the  surface  back  up..  The  near  corner  Is 
somewhat  low,  however  this  corresponds  to  an  area  of  extrapolation.  Figure 
4.3.1.30  shows  a  larger  plot  of  the  surface. 

fjju.re  1.4.4.30.  The  test  surface  Is  shown  In  part  a,  and  part  b  appears 
almost  Indistinguishable  from  It.  There  Is  a  very  slight  flattening  at  the 
right  edge  near  the  center.  Part  c  is  also  a  good  approximation,  with  a  somewhat 
flattened  area  at  the  right,  in  front  of  the  center.  While  there  are  many 
points  r.  urby,  there  Is  a  relatively  large  gap  along  the  edge  which  accounts 
for  poor  definition  of  the  surface  there.  Less  noticeable  Is  a  slightly 
raised  area  to  the  rear  of  the  center,  along  the  right  edge.  The  most  notice¬ 
able  defect  In  part  d  Is  the  poor  behavior  at  the  front  edge,  the  surface 
being  slightly  high  toward  the  right  of  the  center,  and  low  at  the  right  front 


corner. 


Figure  1 .4.5,30.  This  surface  Is  difficult  to  fit  closely  because  of 
Its  sharply  peaked  behavior,  shown  In  part  a.  While  the  peak  In  part  b  is 
well  defined,  to  the  right  one  observes  small  "kinks",  and  generally  wavy 
behavior  around  the  edges,  more  noticeable  at  the  front  and  right  because  of 
the  viewing  point.  The  peak  (at  {•£•,  5-) )  Is  poorly  defined  by  the  set  of 
points  In  part  c  and  thus  1$  considerably  low.  The  wavy  behavior  around  the 
edges  Is  again  observed,  but  somewhat  amplified.  In  part  d  the  peak  Is  higher, 
but  the  wavy  behavior  away  from  the  peak  Is  very  pronounced,  although  the 
surface  Is  smooth  In  the  sense  there  are  no  apparent  "kinks"  In  the  surface. 

Figure  1 .4.6.30.  This  surface,  a  part  of  a  sphere,  Is  shown  In  part  a. 

The  surfaces  shown  In  parts  b,  c,  and  d  show  varying  amounts  of  Imperfect 
behavior,  mostly  appearing  as  flattened  spots  on  the  surface. 

Figure  2.0,0.30.  Part  a  shows  the  cardinal  function.  The  waviness  that 
extends  throughout  the  square  Is  not  desirable  and  Is  probably  an  artifact 
of  the  underlying  polynomial  -  like  Interpolant.  Part  b  shows  the  surface  for 
Aklma's  data  and  basically  appears  reasonable.  Ther?  Is  some  wavy  behavior 
In  the  cross  section  lines  near  the  base  of  the  sharp  rise  toward  the  rear. 

Part  c  shows  a  portion  of  the  surface  for  Ferguson's  data.  Extrapolation  Is  In¬ 
volved  at  the  front  corners  and  the  surface  dips  a  lot  toward  the  right  front. 
The  same  data  scaled  by  a  factor  of  3  In  the  y-directlon  Is  shown  In  part  d. 

Here  the  surface  dips  at  the  left  front  corner  as  well  and  rises  at  the  left 
rear  corner  where  mild  extrapolation  occurs.  Parts  c  and  d  exhibit  some  of 
the  effects  of  scaling  differently  in  the  two  variables. 

Figure  3.1.1.30,  This  figure  shows  the  effects  of  varying  the  parameter 
from  Its  nominal  value.  For  this  method  the  parameter  Is  the  number  of  delta 
Iterates  performed.  The  test  surface  Is  shown  In  part  a.  Part  b  Is  the  surface 
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obtained  with  the  100  point  data  set  after  one  Iteration.  Part  c  shows  the 
surface  after  the  nominal  number  of  Iterations,  which  is  three.  Part  d  shows 
the  surface  after  five  Iterations.  The  surface  shows  definite  Improvement 
with  additional  Iterations,  particularly  when  comparing  b  and  c.  Some 
Improvement  Is  seen  in  d  compared  to  c,  particularly  In  that  the  peak  Is 
filled  out,  although  the  flattened  portion  near  the  left  front  persists. 

The  deviation  statistics  table  P.1  shows  Improvement  In  that  aspect  also. 

Figure  3.3.1.30.  This  figure  Is  the  analogue  of  the  previous,  but  with 
the  25  point  data  set  rather  than  the  100  point  data  set.  In  this  case  It 
Is  not  readily  apparent  that  the  surface  Improves  with  more  Iterates,  although 
part  c  Is  In  some  respects  more  pleasing  than  part  b,  since  the  slight  rise 
toward  the  rear  edge  near  the  rlg^t  corner  has  been  lessened.  .  However  the 
surface  has  been  depressed  at  the  front  edge  near  the  right  corner,  which 
continues  with  more  Iterations,  as  seen  In  part  d.  The  deviation  statistics, 
Table  P.3  show  Improvement  In  the  maximum  deviation,  but  Increases  In  the  mean 
deviation  and  root  mean  square  deviation  after  3  Iterations. 

Figures  4. n. 1.30.  These  figures  are  larger  copies  of  the  figures  given 
in  Figure  1.4.1.30  ,  and  were  previously  discussed. 

Figure  4.1,2,30.  This  function  Is  probably  the  most  difficult  to  fit 
well.  There  Is  some  irregularity  along  the  sharp  rise  about  three-fourths  of 
the  way  back  along  the  diagonal.  The  most  obvious  defects  are  near  the  front 
corner,  on  the  right,  and  a  waviness  along  the  front  edge. 

Figure  4,2.2.30,  The  appearance  of  large  gaps  between  data  points  (see 
Figure  0.2. 0.0)  leads  to  a  quite  wavy  surface.  With  the  exception  of  a  few 
cross  sections  at  about  y  of  the  way  back  along  the  diagonal,  the  surface  Is 
quite  smooth,  however,  with  the  greatest  overshoot  occurlng  near  the  left 


rear  and  right  front  corners,  as  might  be  expected. 

figure  4.3.2.30.  the  regularity  of  the  set  of  data  points  (see 
Figure  O.3.O.O.)  used  here  leads  to  a  more  regular  appearing,  If  somewhat  ' 
wavy,  surface.  The  surface  appears  to  be  vefy  smooth  with  no  apparent  "kinks" 
observable  along  the  cross-sections.  As  usual,  edge  behavior  seems  to  account 
for  the  largest  deviations; 

■Figure  ■4.1 .3. 30.  This  surface  Is  almost  Indistinguishable  from  the 
original.  The  only  apparent  defect  Is  a  slight  flattening  of  the  surface  at 
the  right  edge  near  the  center. 

Figure  4.2.3.30,  The  large  gaps  In  the  data  show  on  this  surface,  but 
less  conspicuously  than  on  some  others.  The  right  end  of  the  surface  appears 
depressed,  with  the  front  edge  also  appearing  to  be  lower  than  the  test  surface. 
The  surface  Is  quite  smooth,  however. 

Figure  4.3.3.30.  This  surface  Is  quite  smooth  and  pleasing,  but  the  left 
rear  corner  Is  considerably  higher  than  that  of  the  test  surface.  The  slope  of 
the  surface  toward  tne  right  center  seems  to  be  more  gentle  than  that  of  the 
test  surface  . 

3.1.0.  Inverse  Distance  Weighted  Methods 

The  performance  of  schemes  within  the  general  class  of  methods  varies  a 
great  deal.  The  basic  Shepard's  method  (program  #18)  with  exponent  2  Is 
unacceptable  for  a  variety  of  reasons  for  all  but  some  very  special  applications 
(perhaps).  For  more  than  a  few  points  the  method  does  not  perform  as  one  would 
be  led  to  believe  when  observing  the  method  for  5  or  10  points.  These  are 
mainly  the  size  examples  given  In  previous  literature.  As  can  be  seen  from  the 
plots,  Figures  1.4.1.18  and  2.0.0.18,  the  surfaces  often  tend  to  have  sharp 


peaks  and  dips  at  the  data  points.  In  fact,  the  resolution  of  the  plots  Is  not 


-66- 


fine  enough  to  show  the  exect  nature  of  the  surface.  . 

Considerable  Improvement  accrues  from  localizing  the  method  (program  #7). 
This  could  likely  be  accomplished  In  a  variety  of  ways  other  than  the  approach 
we  have  taken,  for  example  by  using  an  appropriate  (larger)  exponent.  No 
experimentation  was  done  In  this  direction.  The  plots  for  this  program. 

Figures  1.4. 1.7  and  2. 0.0. 7  show  clearly  the  increased  Influence  of  a  data 
point  on  the  surface  nearby.  This  Is  especially  evident  In  the  cardinal 
function  plot  where  the  Influence  of  a  single  point  Is  seen.  Basically, 
localizing  the  scheme  causes  the  well  advertised  flat  spots  to  become  more 
prominent  and  the  surfaces  to  become  more  pleasing.  However,  for  general, 
purpose  Interpolation  the  scheme  Is  still  basically  unacceptable. 

Forming  boolean  sums  with  other  approximations  does  not  seem  to  work  well 
In  this  case.  It  appears  that  for  the  Idea  of  boolean  sum  approximations  to  work 
the  second  approximation  In  the  boolean  sum  must  be  a  good  approximation  Itself. 
The  least  squares  plane  used  In  our  program  (#2)  Is  not  suitable  since  It  will 
consistently  allow  undershoot  near  peaks  and  hence  appears  to  have  flat  spots 
(not  necessarily  with  zero  slopes)  at  points  where  this  occurs.  This  Is  partic¬ 
ularly  noticeable  In  Figure  1.4. 1.2,  less  so  in  Figure  2. 0.0. 2.  The  alternative 
to  the  least  squares  plane,  a  higher  degree  approximation  such  as  a  quadratic 
would  likely  work  very  well,  but  is  quite  expensive  as  we  will  note  later. 

Another  way  to  generalize  Shepard's  method  Is  through  the  use  of  more 
Information  about  the  surface  near  the  data  points.  Use  of  the  least  squares 
plane  passing  through  each  data  point  In  conjunction  with  a  local  weighting 
function  (program  #3)  leads  to  an  Improved  surface,  however  the  surfaces  often 
tend  to  look  somewhat  lumpy,  as  can  be  seen  from  the  plots.  Figures  1.4. 1.3  and 
2. 0.0. 3.  It  Is  particularly  bad  In  Figure  1.4.1.3c,  probably  because  of  the 


(Intended)  varying  sparseness  of  the  data  points.  The  use  of  planar  fits 
does  not  teem  adequate. 

The  use  of  a  quadratic  function  passing  through  each  data  point  (program 
#17)  Idads  to  VI rtually  no  Improvement  over  the  basic  Shepard's  method  due 
to  the  Influence  of  "far  away"  points.  The  plots  are  shown  In  Figures  1.4.1.17 
and  2. CKO. 17.  and  aire  similar  to  these  for  program  #18.  although  these  are 
generally  somewhat  nicer  In  that  the  number  of  sharp  peaks  Is  reduced. 

The  use  of  a  quadratic  least  squares  fit  at  each  data  point  In  conjunction 
with  localization  of  the  weights  (program  #14)  leads  to  a  significant  Improve¬ 
ment  over  other  methods  of  this  type  In  most  (not  all)  Instances,  especially  for 
larger  numbers  of  data  points.  For  small  numbers  of  data  points  the  surfaces 
seem  to  be  adversely  affected  by  what  might  be  termed  "edge  effects".  This  Is 
not  unique  to  this  scheme  but  occurs  with  other  methods,  particularly  local 
methods.  More  will  be  said  of  this  in  Section  4.  The  plots  shown  In  Figures 
1.4.n.l4,  2.0.0.14,  4.1 .n. 14,  4.2.n.l4,  and  4.3.n.14.,  show  some  of  this,  and 
a  particularly  good  Illustration  Is  given  by  the  cardinal  function  plot, 
2.0.0.14a,  which  can  be  seen  to  behave  In  unseemly  fashion  near  the  left  rear 
corner.  Within  this  class  of  tested  methods,  the  modified  quadratic  Shepard 
method  Is  undoubtedly  the  best  performer  overall.  The  effect  of  changing  the 
parameters  In  the  method  are  shown  In  Figures  3.1.1.14  and  3.3.1.14.  A  greater 
tolerance  to  changes  in  the  parameter  Is  seen  for  th*  larger  data  set,  while 
edge  effects  are  more  prominent  as  the  radius  of  Influence  Is  decreased, 
particularly  for  the  smaller  data  set. 

The  other  approach  to  modifying  Shepard's  method  Is  that  taken  by  McLain 
[39],  and  Implemented  here  In  program  #5  for  the  quadratic  approximation.  From 
the  point  of  view  of  fitting,  the  method  works  quite  well,  although  It  may  be 
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more  prone  to  "edge  effects"  as  can  be  seen  from  the  plots,  Figures  1.4. 1.5 
and  2. 0.0. 5,  especially  2. 0.0. 5  a  and  c.  The  main  reason  for  discounting 
the  method  Is  the  rather  high  computational  burden.  The  modified  quadratic 
Shepard's  method  yields  results  as  good  or  better,  but  at  much  less  cost. 

The  second  McLain  type  of  Interpolant  was  the  fit  with  a  linear  function, 
but  with  modified  weights  (program  #8) .  Because  of  the  necessary  undershoot 
near  peaks  one  can  not  expect  the  method  to  perform  In  a  satisfactory  fashion. 
As  we  see  from  the  plots,  Figures  1.4.1.B  and  2. 0.0. 8,  we  are  not  disappointed 
In  that  respect.  Overall  the  surfaces  appear  to  be  quite  lumpy  and  generally 
unacceptable.  The  cardinal  function  plot  shows  a  somewhat  more  peaked  function 
than  might  be  expected,  almost  like  Shepard's  method. 

3.2.0.  Franke's  Method 

The  performance  of  this  class  of  methods  Is  somewhat  uneven,  giving  quite 
reasonable  results  In  some  cases  and  not  In  others.  Edge  effects  seem  to  come 
Into  the  method  prominently  giving  poor  results  for  data  sets  with  small 
numbers  of  points.  When  the  local  approximations  are  optimal  approximations 
In  (program  #6)  the  plots  are  shown  In  Figures  1.4. 1.6  and  2. 0.0. 6. 

Typical  edge  effects  are  seen  In  Figure  2.0.0.6a,  the  cardinal  function. 

When  the  optimal  approximations  are  taken  boolean  sum  with  a  least 
squares  plane  (program  #1),  the  resulting  surfaces  are  virtually  unchanged, 
as  can  be  seen  In  Figures  1.4. 1.1  and  2. 0.0.1.  Thus  It  appears  that  the  use 
of  the  boolean  sum  with  the  plane  Is  mainly  to  Incorporate  polynomial  pre¬ 
cision,  and  Its  use  with  a  reasonably  good  approximation  will  not  Improve  It 
very  much.  This  observation  was  also  made  by  Foley  In  his  thesis  [15]. 
Additional  plots  are  given  In  Figures  1 . 4 . n . 1 ,  4.1.n.l,  4.2.n.l,  and  4.3.n.l, 
and  basically  show  that  the  method  performs  competently  for  larger  data  sets 
and  not  so  well  on  smaller  ones.  The  variation  of  the  parameter  In  the  metnod, 
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Figures  3. 1.1.1  and  3. 3. 1.1  show  the  effects  of  localizing  too  much  are 


drastic  (plot  b)»  whereas  the  reverse  has  lesser  overall  Impact. 

The  use  of  thin  plate  local  functions  (program  #24)  generally  results 
In  much  more  pleasant  appearing  approximations.  This  Is  particularly  true 
for  smaller  numbers  of  data  points  as  can, be  observed  In  the  plots  ,1n  Figures 

1.4.1.24  and  2.0.0.24.  While  error  statistics  Indicate  some  Improvement 
over  programs  #1  and  #6,  the  plots  generally  appear  to  be  considerably  more 
pleasant.  Additional  surfaces  are  shown  In  Figures  1.4.n.24,  4.1.n.24, 

4.2.n.24,  and  4.3.n.24.  The  varletlon  of  the  parameter  shown  In  Figures 

3.1.1.24  and  3.3.1.24  show  basically  the  same  trend  as  before:  localizing 
too  much  tends  to  degrade  the  approximation,  while  the  reverse  has  less 
Impact.  Overall,  the  use  of  thin  plate  functions  In  the  method  Is  a  nice 
Improvement. 

3.3.0.  Triangle  Based  Blending  Methods 

The  performance  of  this  class  of  methods  Is  dependent  on  the  type  of  nodal 
function  used.  If  they  are  good  local  approximations  to  the  surface,  the 
overall  approximation  will  be  good.  In  line  with  this,  the  linear  nodal 
function  method  (program  #12)  does  not  perform  adequately.  The  plots,  shown 
In  Figures  1.4.1.12  and  2.0.0.12  show  the  transition  between  local  approx¬ 
imations  resulting  In  apparent  creases  In  some  Instances.  This  Is  due  partly 
to  the  resolution  of  the  plots.  Another  defect,  but  one  which  Is  an  artifact 
of  the  triangulation  Is  the  apparent  edge  especially  noticeable  In  Figure 
2.0.0.12c  and  d.  This  Is  due  to  the  occurence  of  a  very  long  thin  triangle 
along  the  edge  of  the  convex  hull.  The  result  Is  that  the  blend  of  approximations 
near  the  middle  part  of  the  triangle  does  not  reflect  the  actual  behavior  at 
nearby  points  (not  in  the  triangle).  This  Is  not  a  defect  unique  to  this 
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method,  but  occurs  to  some  extent  In  all  triangle  based  methods. 

The  use  of  quadratic  nodal  functions  (program  #13)  results  In  a  very 
reasonable  approximation,  as  can  be  seen  In  the  plots  given  In  Figures 
1.4.1.13  and  2.0,0.13.  The  apparent  discontinuity  In  Figure  2.0.0.13d  near 
the  back  1$  due  to  the  previously  mentioned  problem  of  long  thin  triangles 
which  can  occur.  The  other  edge  at  the  left  front  Is  not  nearly  so  apparent 
here.  Additional  plots  are  given  In  Figures  1.4.n.l3,  4.1 .n . 1 3 ,  4.2.n.l3, 
and  4.3. n. 13.  The  parameter  variation  plots  3.1.1.13  and  3.3.1.13  basically 
show  that  too  much  localization  degrades  the  surface.  This  method  generally 
performs  In  quite  an  acceptable  manner  provided  the  disadvantages  of  triangles 
are  acceptable  and  are  outweighed  by  the  advantages  of  triangles  and  the 
overall  method. 

3.4.0.  Finite  Element  Based  Methods 

The  performance  of  a  method  In  this  class  Is  greatly  dependent  on  the 
quality  of  the  estimated  partial  derivatives.  This  Is  the  major  problem  with 
Aklma's  method  (program  #4)  and  causes  the  surfaces  to  have  a  somewhat  lumpy 
and  uneven  appearance.  The  plots  for  Aklma's  method  are  given  In  Figures 
1.4.n.4,  2. 0.0. 4,  4.1.n.4,  4.2.n.4,  and  4.3.n.4.  The  poor  derivative  estimates 
are  especially  noticeable  In  Figure  1.4. 6. 4  where  the  surface  has  a  somewhat 
crumpled  look.  The  variation  of  parameters  plots,  Figures  3. 1.1. 4  and  3. 3. 1.4 
Indicates  that  the  nominal  value  we  have  chosen  Is  probably  about  the  right 
one  to  use.  Figure  3. 3. 1,4  seems  to  show  less  sensitivity  to  the  parameter 
than  3. 1.1 -4.  Figures  2.0.0.4c  and  d  show  the  characteristic  defect  of  triangle 
based  methods.  Aklma's  method  Is  very  fast,  usually  faster  than  other  methods 
by  a  factor  of  4  or  5  or  more. 

Aklma's  method,  modification  one  (program  #10)  performs  somewhat  better 
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than  the  original,,  but, still,  does  not  achieve  the  good  fits  which  the  under- 

\ 

lying  approximations  should  be  capable  of  making.  All  plots,  Figures  1.4.n.l0, 
2.0.0.10,  4 . 1 . n . 1 0 ,  4.2.n.l0,  and  4.3.n.l0  are  basically  similar  to  those  of 
the  original  method,  although  the  statistics  on  deviations  generally  show  it 
to  be  slightly  better. 

Aklma's  method,  modification  two  (program  #11)  would  seem  to  promise  to 
be  better  than  either  of  the  previous,.  However,  the  results  basically  show 
little  or  no  Improvement  over  either  one.  The  plots,  given  In  Figures  1.4.1.11 
and  2.0.0.11  are  basically  similar  to  those  for  programs  #4  and  #10. 

Aklma's  method,  modification  three  (program  #16)  is  undoubtedly  the  best 
performer  of  the  (four.  The  uneveness  of  the  surface  Is  gone  from  nearly  all 
of  the  plots,  shown  In  Figures  1.4.n.l6,  2.0.0.16,  4.1.n.16,  4.2.n.l6,  and 
4.3.n. 16.  Even  the  artifact  of  an  apparent  edge  due  to  the  triangulation, 
usually  prominent  In  the  analogues  to  Figure  2.0.0,16c  and  d,  have  been  reduced 
a  great  deal.  However,  the  cardinal  function  now  shows  some  unbecoming  behavior 
along  the  left  edge  near  the  rear  corner.  This  seems  to  be  caused  by  the 
quadratic  approximation  as  It  also  occurs  In  programs  #13  and  #14,  where  the 
Identical  quadratic  approximation  Is  used.  The  method  seems  to  be  fairly 
insensitive  to  the  parameter  value,  as  Is  shown  In  Figures  3.1.1.16  and  3.3.1.16, 
although  the  larger  value  In  Figure  3 . 1 . 1 . 1 6d  shows  some  more  noticeable 
defects  along  the  front  slope.  Incorporation  of  the  quadratic  to  estimate 
derivatives  results  In  considerably  larger  preprocessing  time. 

Nielson's  minimum  norm  network  (program  #19)  shows  the  capability  of 
triangle  based  approximation,  given  the  appropriate  values  of  the  derivatives. 

As  can  be  seen  In  the  plots  for  the  method,  Figures  1.4.n.l9,  2.0.0.19,  4.1 . n . 1 9 , 
4.2.n.l9,  and  4.3.n.l9,  the  surface  almost  always  appears  quite  smooth  and 


visually  pleasing.  Tho  one  displeasing  behavior  Is  seen  In  Figure  2.0.0.19c 
and  d  (and  to  some  extent  in  b)  where  the  occurence  of  long  slim  triangles 
gives  rise  to  apparent  discontinuities  across  the  triangle.  Again,  this 
simply  reflects  the  fact  that  In  the  middle  of  such  a.  triangle  the  function 
may  not' be  appropriately  represented  since  it  is  far  from  the  vertices  of 
the  triangle,  and  generally  closer  to  other  data  points  which  have  minimal 
Influence.  The  method  Is.  reasonably  fast  and  Is  undoubtedly  the  best  per¬ 
former  In  this  class  of  methods. 

Lawson's  method  (program  #28)  Is  similar  In  spirit  to  Aklma's,  but 
basically  performs  In  much  better  fashion  than  all  but  modification  threr 
of  Aklma's  method  (program  #16).  The  plots,  shown  In  Figures  1 .4. n . 

2.0.0.28,  4.1. n. 28,  4.2.n.28,  and  4.3.n.l6.  One  caution  regarding  the  plots 
for  this  method:  The  program  does  not  extrapolate  outside  the  convex  hull 
of  the  set  {(xk>  yk)l  and  the  function  values  at  such  points  have  been  set 
to  zero  In  the  plots.  Care  should  be  taken  when  viewing  the  plots  based  on 
the  100  point  and  25  point  data  sets,  as  well  as  the  Ferguson  data  set  to  not 
let  these  points  Influence  one's  perception  of  the  surface.  Such  points  are 
omitted  from  the  deviation  statistics  for  this  method.  The  usual  artifact  of 
triangulation  methods  Is  not  seen  in  Figure  2.0.0.28c  and  d  because  no  grid 
points  fall  In  the  appropriate  place.  It  was  extrapolation  In  the  other 
methods  which  made  It  more  visible.  It  Is  noticeable  In  Figure  2,0.0.28b 
along  the  left  edge  and  the  back,  just  as  It  is  in  the  other  triangulation 
based  methods,  Lawson's  method  Is  quite  efficient  in  terms  of  timing,  being 
only  slightly  slower  than  Aklma's  method  (#4). 

3.5.0.  Foley's  Methods 


Foley's  generalized  Newton  Interpol  ant  (program  #25)  performs  in  a  somewhat 


expected  manner  In  that  polynomials  of  high  degree  generally  do  not  work  very  well. 
This  Is  particularly  evident  when  any  extrapolation  Is  Involved,  or  In  regions 
where  no  points  are  nearby.  The  plots  given  In  Figures  1.4.1.25  and  2.0.0.25 
show  this,  especially  In  the  latter,  Figures  2.0.0.25c  and  d,  where  more  extrap¬ 
olation  Is  required  than  In  the  other  examples.  The  cardinal  function,  2.0.0.25a 
has  some  polynomial  type  (mls)behavlor  near  the  front  corner,  also. 

The  Newton  delta  sum  Bernstein  Interpolant  (programs  #26)  Is  something  of 
an  Improvement  In  most  Instances,  although  the  overall  set  of  surface  plots, 
shown  In  Figures  1.4.1.26  and  2.0.0.26  show  basically  the  same  kind  of  behavior. 
The  additional  program  complexity  and  time  required  are  probably  not  worth  the 
result  obtained  here. 

The  use  of  bicubic  splines  and  Iteration  In  connection  with  the  generalized 
Newton  polynomial  (program  #30)  often  results  In  vastly  Improved  surfaces  as 
can  be  seen  from  the  plots  In  Figures  1.4.n.30,  ,2.0.0.30,  4.1 . n . 30 ,  4.2.n.30, 
and  4.3. n. 30.  However,  this  Is  not  universally  true,  and  In  particular  the 
cardinal  function  shown  In  Figure  2.0.0.2a  Is  less  desirable.  Other  surfaces 
based  on  the  25  point  set  are  also  adversely  affected.  This  Is  shown  In  Figure 
3.3.1.30  as  more  Iterations  of  the  delta  sum  produce  poorer  surfaces.  On  the 
other  hand  Figure  3.1.1.30  shows  definite  improvement  as  the  number  of  Iterations 
Increases.  All  things  said,  however,  this  method  seems  to  be  the  best  of  Foley's 
methods.  Computation  time  is  not  excessive,  although  It  Is  slower  than  most 
methods  for  the  problems  In  our  tests. 

The  use  of  a  modified  Shepard's  method  In  place  of  the  generalized  Newton 
Interpolant  In  the  Iterated  delta  sum  with  bicubic  splines  (program  #31)  generally 
gave  poorer  results  than  program  #30.  The  plots  are  shown  in  Figures  1.4.1.31 
and  2.0.0.31.  As  can  be  seen  In  2.0.0.31a,  the  cardinal  function  Is  Improved 


over  that  In  Figure  2.0.0.30a,  as  are  the  surfaces  In  2.0.0.31c  and  d.  That 
the  latter  are  better  Is  no  doubt  due  to  the  stable  extrapolation  of  Shepard's 
method. 

3.6.0.  Global  Basis  Function  Type  Methods 

The  performance  of  this  class  of  methods  varies  widely.  Most  (Duchon's 
methods  are  the  exception)  are  dependent  on  scaling  or  a  parameter  specified 
by  the  user.  We  have  attempted  to  reduce  these  to  an  automatic  value  based  on 
some  estimate  of  mean  distance  between  points. 

The  rotated  Gaussian  (program  #20)  did  not  perform  well.  The  plotSi 
shown  In  Figures  1.4.1.20  and  2.0.0.20  are  quite  smooth  and  give  reasonable 
appearing  approximations  In  the  former.  The  cardinal  function  In  Figure 
2.0.0.20a  appears  to  have  some  undue  Influence  near  the  front  corner.  The 
other  plots  In  Figure  2.0.0.20,  especially  b,  show  a  tendency  of  the  surface 
to  exhibit  local  Guasslan  "bumps".  The  surface  tends  to  zero  far  away  from 
the  data.  Experimentation  with  the  parameter  (related  to  the  variance)  showed 
the  method  to  be  sensitive  to  Its  value,  and  to  depend  on  the  function  values 
rather  than  only  the  (xk,  yk)  points.  For  example,  a  nicer  cardinal  function 
could  be  obtained  by  varying  the  parameter,  but  this  degraded  the  performance 
on  the  surface  shown  In  Figure  1.4.1 . 20d •  which  Is  based  on  the  same  (xk>  yk) 
sets.  For  these  reasons  we  don't  think  this  Is  a  suitable  Idea. 

The  multiquadric  method  proposed  by  Hardy  (program  #21)  performs  very  well. 
The  plots,  shown  In  Figures  1.4,n.21,  2.0.0.21,  4.1 . n . 21 ,  4.2.n.21,  and  4.3.n.21 
show  that  the  method  produces  very  smooth  and  pleasing  surfaces.  The  deviations 
tables.  0,  and  the  related  tables  E  show  that  the  method  Is  consistently  among 
the  most  accurate,  as  well.  The  value  of  the  parameter  was  computed  from  a 
formula  and  Is  related  to  the  mean  distance  to  nearest  neighbor  In  the  set 
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{(x^,  y^)).  The  surface  Is  quite  stable  with  respect  to  changes  In  the  value 
of  the  parameter.  a$  can  be  seen  In  Figures  3.1.1.21  and  3.3,1.21.  The  "best" 
value  Is  probably  dependent  on  the  function  values  as  well,  but  we  obtained 
excellent  results  without  considering  that.  Larger  values  of  the  parameter, 
r,  seemed  to  give  better  results,  but  the  system  of  equations  became  too 
Ill-conditioned  to  solve  In  single  precision.  Thus  a  somewhat  smaller  nominal 
value  was  chosen  than  might  have  been  otherwise.  Larger  values  did  degrade 
the  performance  on  smaller  sets  of  data  while  they  Improved  It  on  larger  sets 
of  data.  The  method  has  also  performed  well  on  recent  tests  with  ocean  bottom 
topography  [47]. 

The  use  of  "reciprocal  multiquadrics"  for  the  basis  functions  (program  #27) 
also  worked  quite  well.  The  plots  are  shown  In  Figures  1.4.n.27,  2.0.0.27, 
4.1.J1.27,  4.2.n.27,  and  4.3.n.27.  The  surfaces  are  again  seen  to  be  very 
smooth.  The  basis  functions  resemble  the  rotated  Guasslan  (#20)  but  generally 
perform  much  more  reliably  than  It.  In  particular,  tne  method  Is  much  less 
sensitive  to  variations  In  the  parameter,  although  very  small  values  of  the 
parameter  will  lead  to  a  surface  consisting  of  sharp  peaks  at  each  data  point 
(or  holes,  If  the  function  value  Is  negative).  For  a  range  of  values  near  the 
nominal  value  chosen  for  the  parameter  the  method  Is  quite  stable.  Overall 
Its  peformance  Is  nearly  as  good  as  the  multiquadric  method. 

3 

The  method  of  Ouchon  which  Involves  the  use  of  basis  functions  dk  (pro¬ 
gram  #22)  works  quite  well.  The  plots,  shown  In  Figures  1.4.1,22  and  2.0.0.22, 
show  very  smooth  surfaces  with  a  pleasing  appearance.  The  cardinal  function. 
Figure  2.0.0.22a  is  very  nicely  shaped.  This  method  was  among  the  better 
performers  overall,  however  solution  of  the  system  of  equations  often  required 
the  use  of  double  precision.  For  this  reason  the  "thin  plate  splines" 
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(program  #23)  were  considered  more  desirable  even  though  the  method  was  not 
superior  In  many  cases.  The  plots  for  the  thin  plate  splines  are  shown  In 
Figures  1.4.n.23,  2.0.0.23,  4.1.n.23,  4.n.2.23,  and  4.n.3.23.  The  surfaces 
are  quite  smooth,  Figures  1.4.1.22  and  1.4.1.23  being  very  similar.  The 
cardinal  function.  Figure  2.0.0.23a,  Is  very  similar  to  that  of  the  previous, 
although  It  seems  to  be  slightly  more  peaked  with  less  undershoot  at  the 
front.  With  the  exception  of  Aklma's  surface.  Figure  2.0.0.23b,  all  calcu¬ 
lations  »  »re  performed  In  single  precision  on  the  IBM  360/67. 

The  use  of  rotated  B-splInes  as  a  basis  function  for  each  data  point 
(program  #29)  yields  variable  results.  The  method  Is  very  sensitive  to  the 
choice  of  radius  at  which  the  function  goes  to  zero.  The  nominal  value 
used  was  chosen  on  the  basis  of  good  performance  for  the  surface  shown  In 
Figure  1.4.1.29.  This  resulted  In  unacceptable  behavior  In  the  cardlr.il 
function.  Figure  2.0.0.29a.  In  this  respect  the  method  seems  to  be  similar 
to  program  #20,  the  rotated  Sausslans.  The  use  of  a  radius  which  resulted 
In  an  acceptable  cardinal  function  seriously  degraded  Its  performance  on  the 
surface  1 1.  Figure  1.4.1.29d,  which  is  based  on  the  same  (x^,  yk)  points.  For 
that  reason  the  method  Is  judged  unacceptable. 

4.0.  Summary 

This  summary  generally  deals  only  with  the  extensively  tested  methods. 

In  the  first  two  groups  of  Table  S.  These  two  groups  were  selected  on  the 
basis  of  meetlngone  or  both  of  two  criteria,  (1)  availability  (are  documented 
programs  readily  available;,,  and  (2)  performance  In  these  tests.  The  following 
local  methods  were  selected  (given  by  number;  refer  to  Table  S  for  a  pointer 
to  the  description):  1,  4,  10,  13,  14,  16,  24,  23.  The  global  methods  selected 
were:  19,  21,  33,  27,  30. 
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The  discussion  of  overall  performance  Is  made  In  separate  sections  for 
local  and  global  methods.  As  we  have  noted  previously,  global  methods  are 
not  feasible  for  very  large  sets  of  data.  If  a  data  set  consists  of  some 
100  -  200  points.  It  Is  feasible  and  practical  to  use  a  global  method.  Me 
have  demonstrated  here  that  for  100  points  or  less,  some  real  advantages 
accrue  for  global  methods.  In  particular,  global  methods  perform  better  In 
terms  of  their  deviations  from  test  surfaces.  This  Is  seen  In  Table  E.2, 
where  local  methods  are  "best"  In  only  4  of  18  cases.  Global  methods  seem 
to  be  less  likely  to  exhibit  edge  effects  than  local  methods.  It  Is  possible 
that  this  Is  partially  responsible  for  the  results  In  Table  E.2,  although  the 
surfaces  for  global  methods  are  generally  much  smoother  and  more  pleasant 
appearing  than  those  for  local  methods.  For  very  large  sets  of  data  the 
regions  where  edge  effects  occur  should  be  a  smaller  part  of  the  overall 
region  of  Interest. 

The  processing  time  for  local  methods  Is  generally  less  than  that  for 
global  methods,  although  some  global  methods  are  faster  than  some  local 
methods.  The  trend  of  local  methods  Is  generally  to  Increase  at  least  linearly 
with  the  number  of  points  when  total  time  Is  considered.  Global  methods  gener¬ 
ally  Increase  at  least  linearly  also  since  all  points  must  be  Inspected  The 
potential  savings  for  local  methods  comes  from  not  having  to  solve  a  linear 
system  of  N  equations  In  the  preprocessing  phase  and  use  of  only  nearby  points 
In  the  evaluation  phase.  Specific  comments  are  made  In  the  next  two  sections. 
4.1 .  Local  Methods 

The  best  performing  local  methods  are  probably  the  Modified  Quadratic 
Shepard  Method  and  a  similar  program  based  on  a  triangulation  of  the  convex 
hull,  the  Nlelson-Franke  Quadratic  Triangle  method.  Both  methods  perform 
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consistently  well,  with  the  method  based  on  triangles  generally  being  slightly 
more  accurate  and  faster,  it  Is,  however  somewhat  prone  to  misbehave  when 
long  slim  triangles  occur,  and  the  triangulation  requires  a  great  deal  of 
auxiliary  storage.  For  general  purpose  use  the  Modified  Quadratic  Shepard's 
Method  Is  favored  for  several  reasons.  It  Is  easy  to  Implement  and  generalizes' 

In  rather  easy  fashion  to  higher  dimensional  spaces.  The  apparent  time 
penalty  In  the  evaluation  phase  could  be  reduced  by  some  additional  preprocessing 
and  auxiliary  storage  to  allow  quicker  determination  of  points  which  (potentially) 
affect  the  Interpolant.  This  sort  of  scheme  should  be  Incorporated  when  dealing 
with  very  large  sets  of  data  to  avoid  excessive  evaluation  times.  Triangle 
based  methods  (whether  local  or  global)  have  this  kind  of  scheme  bull*  In  and 
partly  accounts  for  their  efficiency  during  the  evaluation  phase. 

Aklma's  method  suffers  from  poor  estimates  of  the  derivative  values. 
Modification  one  results  In  modest  Improvement,  generally,  but  still  does  not 
perform  as  well  as  is  possible  for  the  underlying  approximation.  Both  versions 
are  subject  to  the  appearance  of  extraneous  bumps,  as  seen  In  Figures  4. 1.6. 4 
and  4.1.6.10,  as  well  as  overshoot  as  seen  In  Figures  1.4.1.4b  and  1.4.1.10b. 
Lawson's  method,  based  on  a  similar  Idea,  but  with  a  different  element  and 
tetter  estimates  of  the  derivatives,  generally  performs  better  than  Aklma's 
method.  The  use  of  Inverse  distance  weighted  least  squares  quadratics  to 
estimate  the  derivatives  in  Aklma's  method  generally  (but  not  always)  results 
In  Improved  accuracy,  and  even  more  often  gives  a  much  more  pleasing  surface, 
as  can  be  seen  by  comparing  the  corresponding  plots  for  the  three  different 
schemes  for  estimating  derivatives.  This  latter  version  of  Aklma's  program 
performs  about  on  a  par  with  Lawson's  method.  It  requires  considerably  more 
preprocessing  time,  however.  It  Is  subject  to  edge  effects  In  some  cases. 
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Lawson ' s  method  generally,  gives  smooth  appearing  surfaces,  although 
the  occuranqe  of  long  slim  triangles  can  cause  problems.  Some  of  these 
effects  are  not  apparent  because  no  evaluation  points  fall  within  the  long 
slim  triangles  that  occur  along  edges,  as  In  Ferguson's  data,  for  example. 

While  1£  Is  not  a  local  method,  It  should  be  mentioned  that  Nielson's 
Minimum  Norm  Network, method  generally  performs  better  than  any  of  the  local 
methods  based  on  triangles.  It  avoids  the  usual  storage  problem  Involved  in 
solving  a  large  system  of  equations  by  solving  Iteratively,  and  rapid  conver¬ 
gence  Is  obtained.  The  overall  storage  requirements  are  similar  to  Aklma's 
method  (and  Its  variants)  but  more  than  for  Lawson's  method.  Timing  Is  some¬ 
what  slower  than  all  but  modification  three  to  Aklma's  method  In  the  preprocessing 
phase.  It  Is  also  slower  In  the  evaluation  phase,  but  a  different  Implementation 
of  this  could  result  In  It  being  about  as  fast  as  Lawson's  method.  The  under¬ 
lying  approximation  Is  somewhat  more  complicated  than  Aklma's,  but  use  of  an 
evaluation  phase  following  a  strategy  similar  to  Aklma's  should  not  be  slower 
by  a  factor  of  more  than  about  two.  So  far  this  modification  has  not  been  made. 

The  remaining  two  local  methods  are  due  to  the  Investigator.  The  under¬ 
lying  Idea  of  partitioning  the  plane  Into  rectangles  seems  to  be  sound,  result¬ 
ing  In  reasonable  (not  fast,  but  nearly  Independent  of  N)  evaluation  times. 

The  use  of  thin  plate  splines  as  the  local  approximations  Is  a  definite  Improve¬ 
ment  In  both  the  appearance  and  accuracy  of  the  method.  Overall,  however, 
performance  of  the  method  Is  somewhat  disappointing.  It  seems  there  Is  no 
Inherent  reason  why  Its  performance  should  not  be  nearly  as  good  as  the  local 
approximations,  which  are  very  good,  according  to  our  results  on  global  methods. 

It  seems  that  the  amount  of  overlap  In  the  local  Interpolating  functions  Is  not 
sufficient  to  prevent  transition  from  one  rectangle  to  another  resulting  in 
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transition  to  a  fundamentally  different  surface,  as  can  be  seen,  for 
W  example,  In  Figure  1.4.5.24c.  Use  of  much  larger  values  of  NPPR  could 

result  In  better  approximations,  although  our  tests  (In  Tables  P.1  and  P.3.) 

% 

show  conflicting  evidence.  No  further  experimentation  has  been  performed.' 

The  amount  of  auxiliary  storage  required  Is  mild,  compared  to  the  methods  based 
on  triangles,  particularly.  Evaluation  times  are  nearly  Independent  of  the 
number  of  data  points,  which  could  be  useful  for  very  large  N.  Triangle 
based  methods  also  possess  this  property  and  are  faster  than  rectangle  based 
methods  because  faster  evaluation  of  the  local  Interpolants  Is  possible. 

All  things  considered,  the  method  of  choice  here  seems  to  be  the 
Modified  Quadratic  Shepard's  Method.  Its  advantages  of  simplicity  and  mild 
auxiliary  storage  requirements  overcome  Its  relatively  expensive  evaluation 

|  phase.  A  preprocessing  phase  to  determine  (potential)  points  which  affect 

, 

the  Interpolant  In  various  regions  could  be  Implemented  at  a  modest  cost  In 
^  '  time  (probably  less  than  one  second)  and  storage  (about  5/3  N  locations), 

which  for  large  data  sets  would  probably  result  in  evaluation  times  of  10 
seconds  or  less  (Independent  of  N),  but  this  has  not  been  implemented  as  yet. 
4.2.  Global  Methods 

The  most  Impressive  method  In  these  tests  is  the  multiquadric  method  of 
Hardy.  It  Is  consistently  best  or  near  best  Iri  terms  of  accuracy,  and  always 
results  In  visually  pleasant  surfaces.  Nonetheless,  a  certain  skepticism 
persists  because  the  method  has  no  apparent  mathematical  basis  to  explain  Its 
efficacy.  In  some  respects  the  basis  functions  are  somewhat  similar  to  the 
thin  plate  splines  of  Duchon  In  that  they  take  on  large  values  at  points  far 
away  from  the  data  point.  Further,  they  appear  (for  r  »  0)  to  fit  the  class 
of  approximations  discussed  by  Melnguet  [41],  but  the  proofs  do  not  hold. 
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In  the  degenerate  case  (r  ■  0),  Initially  Investigated  by  Hardy,  the  multi- 
quadrics  are  cdnes  with  zero  value  at  the  data  point  (just  as  d  log  d  Is  zero 
at  the  data  point).  Given  only  the  basis  functions  d  log  d  one  might  also  be 
perplexed  at  how  well  they  work.  Perhaps  there  Is  an  equally  elegant  (but 
unknown)  theory  to  explain  the  abilities  of  the  multlquadrlc  method.  On  the 
basis  of  our  tests  we  can  recommend  the  use  of  either  the  multlquadrlc  method 
or  the  thin  plane  splines  as  the  best  of  the  global  basis  function  methods, 
and  perhaps  the  best  of  all  global  methods  considered.  The  reciprocal  multl¬ 
quadrlc  has  some  potential  bad  effects  and  for  too  small  a  value  of  r  will 
give  poor  results,  as  noted  In  Section  3.6.0.  There  seems  to  be  no  reason  to 
use  It  rather  than  the  multlquadrlc  method. 

Nielson's  minimum  norm  network  has  been  discussed  a  bit  In  connection 
with  local  methods  In  the  previous  section.  Computationally  It  Is  a  viable 
method  for  larger  sets  of  data  than  the  methods  requiring  solution  of  a  full 
system  of  N  or  more  equations  since  It  uses  Iteration  on  a  sparse  system  of 
equations.  It  does  use  considerable’ storage  which  will  probably  limit  the 
method  before  excessive  computation  time.  The  use  of  the  method  must  be  done 
with  the  knowledge  that  poor  behavior  can  occur  In  long  slim  triangles,  a 
caution  that  applies  to  all  methods  based  on  triangles.  Nielson's  method  is 
reportedly  soon  to  be  available  In  a  version  which  does  not  require  a  convex 
region,  and  this  could  easily  be  used  to  eliminate  undesirable  triangles 
along  the  edges.  Extrapolation  will  not  be  so  easily  Implemented  In  this 
version,  however,  so  If  that  Is  Important,  It  Is  a  consideration. 

Foley's  TF  delta  sum  bicubic  spline  Is  a  relatively  poor  performer  here. 
Results  using  the  method  have  been  discussed  in  great  detail  as  an  example  In 
Section  3.0.0.  While  the  method  yields  some  very  nice  Interpolants,  It  Is 


-82- 


rather  inconsistent  and  often  has  undesirable  ripples. 

5.0.  Epilog 

This  Investigation  has  consumed  a  great  deal  of  time  and  effort. 

Thanks  are  due  to  numerous  colleagues,  among  them  Greg  Nielson.  Bob  Barnhill, 
Frank  Little,  Tom  Foley,  Rosemary  Chang,  and  others  with  whom  I  discussed 
many  ideas  and  who  made  valuable  (sometimes  followed!)  suggestions.  Thanks 
are  also  due  to  those  who  supplied  working  programs,  among  them  Greg  Nielson, 
Hiroshi  Aklma,  Charles  Lawson,  and  Tom  Foley.  Last,  but  hardly  least,  thanks 
go  to  Linda  Dent  for  her  patience  and  good  humor  during  the  period  when  the 
manuscript  was  being  typed  and  revised,  and  (especially)  during  the  process  of 
pasting  up  plots, 

Despite  the  number  of  Ideas  explored  and  programs  written  and  tested, 
there  are  a  number  of  Ideas  which  were  not  Investigated.  Among  them  are 
several  from  the  CAGD  group  at  the  University  of  Utah.  Many  of  these  are 
based  on  the  use  of  trl angulations,  which  the  Investigator  feels  are  much 
more  suitable  for  the  design  problem  (where  long  slim  triangles  can  be 
avoided)  than  for  the  scattered  data  Interpolation  problem.  It  was  not 
possible  to  test  Vlttltow's  [55]  variations  of  Maude's  method  [37],  although 
It  appears  they  may  perform  reasonably  well.  Another  Idea  which  was  not  tested 
has  its  genesis  In  Briggs  [8],  and  is  available  commercially  [59],  The  user's 
manual  contains  some  Impressive  material,  but  no  tests  of  the  software  have 
been  conducted.  There  are  no  doubt  more  Ideas  worthy  of  investigation  appearing 
In  the  literature. 

In  terms  of  the  data  considered  here.  It  was  for  the  most  part  rather 
nice  data,  even  though  some  effort  was  made  to  Include  some  data  with  varying 
densities.  Real  data  exists  which  Is  very  sparce  In  certain  regions,  or 


lies  In  clumps.  Some  methods  will  not  work  In  a  reasonable  fashion  for  this 
type  of  data,  although  we  have  not  tried  to  determine  which  methods  will  and 
which  will  not.  Methods  based  on  quadratic  approximations  will  likely  mis¬ 
behave  for  such  data.  In  addition,  local  methods  based  on  distance  weighting 
may  have  holes  In  the  domain  of  definition  when  density  varies  greatly  or 
when  data  appears  in  clumps.  Some  additional  work  Is  necessary  to  see  If 
there  are  suitable  local  methods  for  such  data. 

The  Investigator  Is  willing  to  make  further  tests  (at  least  for  the 
supplier  and  perhaps  for  wider  dissemination)  of  working  programs,  under  the 
following  (negotiable)  guidelines:  (1)  The  program  Is  to  be  supplied  on 
cards  (preferably  EBCDIC  punch).  (2)  The  program  Is  to  be  In  the  form  of 
one  or  more  subroutines,  and  a  grid  of  Interpol ant  values  Is  to  be  returned 
by  calling  one  of  them  with  the  appropriate  data  and  workspace.  (3)  The 
program  Is  to  be  In  ANSI  standard  Fortran.  (4)  Program  documentation  will 
be  supplied  to  enable  use  of  the  program.  (S)  A  sample  driver  program  will  be 
supplied.  The  Investigator  will  supply  at  least  the  plots  of  the  type 
1.4.1,n  and  2.0.0.n  and  the  corresponding  error  and  timing  statistics. 

Depending  on  those  results,  additional  tests  may  be  performed  and  reported 
to  the  supplier. 
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Table  5 
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(■ 

I 

l 
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Method 

Maximum 

Mean 

RMS 

Deviation 

Deviation 

Deviation 

1 

Franke  -  3 

.0919 

.00842 

.0148 

'Vjj 

4 

Aklma 

.0647 

.00787 

.0125 

10 

Aklma  Mod.  I 

.0856 

.00784 

.0133 

yi 

13 

Nielson  -  Franke  Q 

.0782 

.00741 

.0122 

3 

14 

Mod.  Quad.  Shepard 

.0573 

.00785 

.0128 

i 

16 

Aklma  Mod.  Ill 

.0520 

.00729 

.0117 

■fi 

24 

Franke  -  TPS 

.0940 

.00887 

.0164 

""if 

28 

Lawson 

.0951 

.00783 

.0124 

19:  Nielson  MlnNorm 

.0492 

.00537 

.00940 

'V 

21 

Hardy  Quadric 

.0225 

.00181 

.00357 

23 

Duchon  TPS 

.0518 

.00525 

.00947 

27 

Hardy  Reclp.  Quad. 

.0247 

.00283 

.00518 

30:  Foley  III 

.0636 

.00473 

.00941 

2 

Mod.  Shepard  ®  Plane 

.156 

.0137 

.0254 

i 

3 

Mod.  Linear  Shepard 

.104 

.00982 

.0172 

■ 

5 

McLain  Mir, 

Franke -l10 

.0601 

.00747 

.0124 

6 

.108 

.0103 

.0188 

i 

7 

Mod.  Shepard 

.224 

.0272 

.0440 

8 

Mod.  McLain  Mg 

.194 

.0167 

.0316 

11 

Aklma  Mod.  II 

.105 

.00875 

.0152 

0 

12 

Nielson  -  Franke  L. 

.125 

.0101 

.0189 

i 

17 

Quad  Shepard 

.264 

.0396 

.0594 

18 

Shepard 

.273 

.0417 

.0620 

20 

Rotated  Gausslans 

.0624 

.00599 

.0112 

22 

Duchon 

.0247 

.00311 

.00578 

25 

Foley  I 

.201 

.0153 

.0305 

26 

Foley  II 

.144 

.0120 

.0229 

29 

Rotated  B-Spllnes 

.0488 

.00790 

.0112 

: 

31 

Foley  IV 

.128 

.0113 

.0204 

Deviations  from  Exponential  test  surface,  100  points 

Table  D.1.1 


-*^D,«r^r  n?  *****  '►to*4-**  v 


•<ut  n»« rtuf*‘ 


“Jf-r: 


Method 


1:  Franke  -  3 
4:  Aklma 
10:  Aklma  Mod.  1 
13:  Nielson  -  Franke  Q 
14:  Mod.  Quad.  Shepard 
16:  Aklma  Mod.  Ill 
24:  Franke  -  TPS 
28 :  Lawson 


19: 

21: 

23: 

27: 

30: 


Nielson  MlnNorm 
Hardy  quadric 
Duchon  -  TPS 
Hardy  Reclp.  Quad. 
Foley  III 


Maximum 

Deviation 

.0518 

.0520 

.0473 

.0721 

.0468 

.0958 

.0295 

.0280 

.0424 

.0244 

.0344 

.0379 

.0281 


Mean 

Deviation 

.00286 

.00303 

.00257 

.00268 

.00264 

.00293 

.00243 

.00221 

.00181 

.00177 

.00210 

.00192 

.00223 


RMS 

Deviation 

.00586 

.00609 

.00542 

.00683 

.00551 

.00809 

.00483 

.00448 

.00434 

.00330 

.00436 

.00388 

.00419 


Deviations  from  Cliff  test  surface,  100  points 
Table  D.1.2 


Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1 :  Franke  -  3 

.0198 

.00164 

.00294 

4:  Aklma 

.0274 

.00224 

.00423 

10:  Aklma  Mod.  1 

.0254 

.00198 

.00367 

13:  Nielson  -  Franke  Q 

.0168 

.00110 

.00206 

14:  Mod.  Quad.  Shepard 

.0125 

.00112 

.00194 

16:  Aklma  Mod.  Ill 

.0142 

.00105 

.00202 

24:  Franke  -  TPS 

.0165 

.00157 

.00273 

28:  Lawson 

.0565 

.00149 

.00359 

19:  Nielson  MlnNorm 

.0195 

.00091 

.00200 

21:  Hardy  Quadric 

.00461 

.00025 

.00052 

23:  Duchon  -  TPS 

.00597 

.00049 

.00092 

27:  Hardy  Reclp.  Quad. 

.00928 

.00068 

.00136 

30:  Foley  III 

.0117 

.00117 

.00196 

Deviations  from  Saddle  test  surface,  100  points 


Method 


Maximum 

Deviation 


Mean 

Deviation 


RMS 

Deviation 


1 :  Franke  -  3 
4:  Aklma 
10:  Aklma  Mod.  I 
13:  Nielson  -  Franke  Q 
14:  Mod.  Quad.  Shepard 
16:  Aklma  Mod.  Ill 
24:  Franke  -  TPS 
28:  Lawson 

19:  Nielson  MinNorm 
21 :  Hardy  Quadric 
23:  Duchon  -  TPS 
27:  Hardy  Reclp.  Quad. 
30:  Foley  III 


.0114 

.0101 

.00675 

.00517 

.00388 

.00330 

.00560 

.00899 

.00303 

.00102 

.00294 

.00227 

.00604 


.00122 

.00124 

.00102 

.00058 

.00065 

.00049 

.00103 

.00061 

.00047 

.00005 

.00017 

.00034 

.00083 


.00189 

.00177 

.00143 

.00083 

.00089 

.00070 

.00141 

.00109 

.00069 

.00011 

.00030 

.00050 

.00117 


Deviations  from  Gentle  test  surface,  100  points 


3  -  ,  -■ 


. 

1 

r 

f. 

? 

Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1; 

1: 

Franke  -  3 

.0358 

.00228 

.00447 

}' 

4: 

Aklma 

.0434 

.00242 

.00510 

r 

1 

10: 

Aklma  Mod.  I 

.0317 

.00215 

.00436 

!' 

13: 

Nielson  -  Franke  0 

.0206 

.00176 

.00337 

l 

14: 

Mod.  Quad.  Shepard 

.0218 

.00182 

.00361 

t 

16: 

Aklma  Mod.  Ill 

.0212 

.00171 

.00337 

1 

24: 

Franke  -  TPS 

.0284 

.00212 

.00418 

E 

{ 

28: 

Lawson 

.0216 

.00154 

.00323 

<L 

ii 

19: 

Nielson  MlnNorm 

.0195 

.00101 

.00229 

J 

21: 

Hardy  Quadric 

.00280 

.00012 

.00031 

{  ( 

f- 

23: 

Duchon  -  TPS 

.0175 

.00088 

.00217 

27: 

Hardy  Reclp.  Quad. 

.00736 

.00030 

.00078 

1 

30: 

Foley  III 

.0143 

.00172 

.00282 

Deviations  from  Steep  test  surface,  100  points 
Table  D.l .5 
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Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1: 

Franke  -  3 

.0119 

.00126 

.00206 

4: 

Aklma 

.0196 

.00196 

.00313 

10: 

Aklma  Mod.  I 

.0172 

.00173' 

.00286 

13: 

Nielson  -  Franke  Q 

.00343 

.00022 

.00043 

14: 

Mod.  Quad.  Shepard 

.00361 

.00026 

.00050 

16: 

Aklma  Mod.  Ill 

.00796 

.00058 

.00094 

24: 

Franke  -  TPS 

.0111 

.00138 

.00206 

28; 

Lawson 

.00954 

.00038 

.00099 

19: 

Nielson  MlnNorm 

.0117 

.00077 

.00165 

21: 

Hardy  Quadric 

.0106 

.00041 

.00111 

23: 

Duchon  -  TPS 

.0170 

.00053 

.00150 

27: 

Hardy  Reclp.  Quad. 

.0241 

.00117 

.00263 

30: 

Foley  III 

.00965 

.00127 

.00203 

Deviations  from  Sphere  test  surface,  100  points 
Table  D.1.6 
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Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1 

Franke  -  3 

'  .347 

.0477 

.0732 

4 

Aklma 

.158 

.0384 

.0535 

10 

Aklma  Mod.  I 

.197 

.0400 

.0570 

13 

Nielson  -  Franke  Q 

.150 

.0326 

.0455 

14 

Mod.  Quad.  Shepard 

.184 

.0340 

.0478 

16 

Aklma  Mod.  Ill 

.164 

.0372 

.0521 

24 

Franke  -  TPS 

.218 

.0346 

.0517 

28 

Lawson 

.287 

.0462 

.0657 

19 

Nielson  MinNorm 

.150 

.0305 

.0437 

21 

Hardy  Quadric 

.137 

.0181 

.0269 

23 

Duchon  -  TPS 

.153 

.0293 

.0421 

27 

Hardy  Reclp.  Quad. 

.140 

.0153 

.0244 

30 

Foley  III 

"  .296 

.0350 

.0546 

2 

Mod.  Shepard  9  Plane 

.208 

.0402 

.0604 

3 

Mod.  Linear  Shepard 

.321 

.0566 

.0870 

5 

McLain  M1Q 

.217 

.0438 

.0625 

6 

Franke  -  l 

"  .357 

.0484 

.0741 

7 

Mod.  Shepard 

.377 

.0571 

.0872 

8 

Mod.  McLain  Mg 

V  \193 

.0379 

.0566 

11 

Aklma  Mod.  II  , 

,  .232 

.0401 

.0582 

12 

Nielson  -  Franke  L. 

.274 

.0446 

.0651 

17 

Quad.  Shepard 

.223 

.0701 

.0915 

18 

Shepard 

.225 

.0709 

.0922 

20 

Rotated  Gausslans 

.137 

.0174 

.0287 

22 

Duchon 

.140 

.0235 

.0338 

25 

Foley  I 

.162 

.0277 

.0387 

26 

Foley  II 

.161 

.0281 

.0383 

29 

Rotated  B-Spllnes 

.137 

.0210 

.0337 

31 

Foley  IV 

.273 

.0422 

.0626 

Deviations  from  Exponential  Test  Surface,  33  points 


c 


Method 

Maximum 

Mean 

RMS 

Deviation 

Deviation 

Deviation 

1: 

Franke  -  3 

.0776 

.0124 

.0190 

4: 

Aklma 

.0543 

.00850 

.0133 

10: 

Aklma  Mod.  I 

.0518 

.00747 

.0122 

13: 

Nielson  -  Franke  Q 

.0878 

.0137 

.0219 

14: 

Mod.  Quad.  Shepard 

.0876 

.0121 

.0206 

16: 

Aklma  Mod.  Ill 

.0680 

.0106 

.0176 

24: 

Franke  -  TPS 

.0561 

.00913 

.0147 

28: 

Lawson 

.0956 

.0126 

.0205 

19: 

Nielson  MlnNorm 

.0582 

.00800 

.0140 

21: 

Hardy  Quadric 

.0577 

.0129 

.0170 

23: 

Duchon  -  TPS 

.0526 

.00777 

.0134 

27: 

Hardy  Reclp.  Quad. 

.0500 

.00853 

.0130 

30: 

Foley  III 

.0914 

.0165 

.0262 

Deviations  from  Cliff  test  surface,  33  points 


Table  D.2.2 


Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1 :  Franke  -  3 

.111 

.0121 

.0224 

4:  Aklma 

.0578 

.0110 

.0165 

10:  Aklma  Mod.  I 

.0578 

.0104 

.0156 

13:  Nielson  -  Franke  Q 

.0679 

.00939 

.0146 

14:  Mod,  Quad.  Shepard 

.0724 

.00907 

.0139 

16:  Aklma  Mod.  Ill 

.0597 

.0104 

.0162 

24:  Franke  -  TPS 

.0662 

.0109 

.0175 

28:  Lawson 

.0685 

.0133 

.0199 

19:  Nielson  MlnNorm 

.0571 

.0102 

.0159 

21 :  Hardy  Quadric 

.0262 

.00442 

.00689 

23:  Duchon  -  TPS 

.0574 

.00912 

.0140 

27:  Hardy  Reclp.  Quad. 

.0505 

.00571 

.00970 

30:  Foley  III 

.0885 

.00888 

.0148 

Deviations  from  Saddle  test  surface,  33  points 


Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1:  Franke  -  3 

.0446 

.00608 

.0101 

4:  Aklma 

.0167 

.00487 

.00623 

10:  Aklma  Mod.  I 

.0160 

.00442 

.00573 

13:  Nielson  -  Franke  Q 

.0312 

.00422 

.00637 

14:  Mod.  quad.  Shepard 

.0272 

.00451 

.00679 

16:  Aklma  Mod.  Ill 

.0204 

.00394 

.00565 

24:  Franke  -  TPS 

.0339 

.00681 

.0107 

28:  Lawson 

.0269 

.00552 

.00815 

19:  Nielson  MlnNorm 

.0214 

.00371 

.00563 

21:  Hardy  Quadric 

.00724 

.00121 

.00204 

23:  Duchon  -  TPS 

.0259 

.00415 

.00714 

27:  Hardy  Reclp.  Quad. 

.0188 

.00266 

.00485 

30:  Foley  III 

.0349 

.00438 

.00674 

Deviations  from  Gentle  test  surface,  33  points 
Table  D.2.4 
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Method 

Maximum 

Deviation 

Mean 

Deviation 

RMS 

Deviation 

1: 

Franke  -  3 

.143 

.0162 

.0298 

4: 

Aklma 

.115 

.0120 

.0240 

10: 

Aklma  Mod.  I 

.109 

.0113 

.0227 

13: 

Nielson  -  Franke  Q 

.0835 

.0104 

.0181 

14: 

Mod.  quad.  Shepard 

.110 

.0113 

.0220 

16: 

Aklma  Mod.  Ill 

.115 

.0119 

.0240 

24: 

Franke  tps 

.150 

.0148 

.0305 

28: 

Lawson 

.139 

.0129 

.0289 

19: 

Nielson  MlnNorm 

.115 

.0106 

.0228 

21: 

Hardy  Quadric 

.0716 

.00850 

.0148 

23; 

Duchon  -  TPS 

.149 

.01 30 

.0296 

27: 

Hardy  Reclp.  quad 

.0963 

.00878 

.0180 

30: 

Foley  III 

.110 

.0143 

.0249 

Deviations  from  Steep  test  surface,  33  points 


Table  D.2,5 


Method 


Maximum 

Deviation 


Mean 

Deviation 


RMS 

Deviation 


12 

Franke  -  3 

.0464 

.00614 

.0106 

4: 

Akima 

.0383 

.00796 

.0110 

10: 

Akima  Mod.  I 

.0393 

.00732 

.0104 

13: 

Nielson  -  Franke  Q. 

.0983 

.00585 

.0177 

14: 

Mod.  Quad.  Shepard 

.101 

.00400 

.0136 

16: 

Akima  Mod.  Ill 

.0819 

.00556 

.0139 

24: 

Franke  -  TPS 

.0307 

.00629 

.00886 

28: 

Lawson 

.0137 

.00210 

.00313 

19: 

Nielson  MlnNorm 

.0186 

.00273 

.00460 

21: 

Hardy  Quadric 

.0203 

.00278 

.00473 

23: 

Duchon  -  TPS 

.0232 

.00315 

.00545 

27: 

Hardy  Reclp.  Quad. 

.0351 

.00414 

.00737 

30: 

Foley  III 

.0269 

.00493 

.00726 

Deviations  from  Sphere  test  surface,  33  points 


Table  D.2.6 
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Method 

Maximum 

Mean 

Deviation 

Deviation 

1 

Franke  -  3 

.240 

.0359 

4 

Aklma 

.134 

.0282 

10 

Aklma  Mod.  I 

.129 

.0280 

13 

Nielson  -  Franke  Q 

.153 

.0350 

14 

Mod.  Quad.  Shepard 

.158 

.0353 

16 

Aklma  Mod.  Ill 

.155 

.0355 

24 

Franke  -  TPS 

.129 

.0267 

28:  Lawson 

.202 

.0327 

19 

Nielson  MlnNorm 

.124 

.0235 

21 

Hardy  Quadric 

.119 

.0235 

23 

Duchdn  TPS 

.121 

.0253 

27 

Hardy  Recip.  Quad. 

.119 

.0214 

30 

Foley  III 

.165 

.0196 

2 

Mod.  Shepard  #  Plane 

.167 

.0328 

3 

Mod.  Linear  Shepard 

.254 

.0418 

5 

McLain  M10 

.255 

.0369 

6 

Franke  -  1 

.241 

.0356 

7 

Mod.  Shepard 

.212 

.0481 

8 

Mod.  McLain  Mg 

.262 

.0377 

11 

Aklma  Mod.  II 

.126 

.0284 

12 

Nielson  -  Franke  L. 

.249 

.0366 

17 

Quad.  Shepard 

.233 

.0550 

18 

Shepard 

.238 

.0559 

20 

Rotated  Gausslans 

.118 

.0237 

22 

Duchon 

.117 

.0246 

25 

Foley  I 

.200 

.0375 

26 

Foley  II 

.166 

.0333 

29 

Rotated  B-Spllnes 

.131 

.0279 

31 

Foley  IV 

.121 

.0195 

Deviations  from  Exponential  Test  Surface, 

25  points 

Table  D.3.1 


C 


RMS 

Deviation 

.0486 

.0386 

.0390 

.0478 

.0486 

.0484 

.0374 

.0458 

.0328 

.0322 

.0348 

.0294 

.0310 


.0466 

.0593 

.0529 

.0484 

.0661 

.0579 

.0396 

.0513 

.0670 

.0709 

.0321 

.0330 

.0517 

.0449 

.0368 

.0276 
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Method 

Maximum 

Mean 

RMS 

Deviation 

Deviation 

Deviation 

Is 

Franke  -  3 

.161 

.0225 

.0408 

4: 

Akima 

.0999 

.0148 

.0257 

10: 

Aklma  Mod.  I 

.0987 

.0143 

.0252 

13: 

Nielson  -  Franke  Q 

.148 

.0166 

.0304 

14: 

Mod.  Quad.  Shepard 

.163 

.0166 

.0314 

16: 

Aklma  Mod.  Ill 

.146 

.0164 

.0305 

24: 

Franke  -  TPS 

.106 

.0148 

.0257 

2B: 

Lawson 

.132 

.0164 

.0283 

19: 

Nielson  MlnNorm 

.0942 

.0138 

.0242 

21: 

Hardy  Quadric 

.0995 

.0143 

.0231 

23: 

Duchon  -  TPS 

.101 

.0135 

.0235 

27: 

Hardy  Reclp.  Quad 

.105 

.0139 

.0236 

30: 

Foley  III 

.0832 

.0165 

.0250 

Deviations 

from  Cliff  test  surface, 

25  points 

Table  D.3.2 


c 

-107- 


Method 

Maximum 

Deviation 

Mean 

Deviation 

1: 

Franke  -  3 

.0688 

.0111 

4: 

Aklma 

.0864 

.0121 

10: 

Aklma  Mod.  I 

.0866 

.0119 

13: 

Nielson  -  Franke  Q 

.0794 

.0115 

14: 

Mod.  Quad.  Shepard 

.0759 

.0114 

16: 

Aklma  Mod.  Ill 

.0787 

.0116 

24: 

Franke  -  TPS 

.0714 

.00983 

28: 

Lawson 

.0875 

.0126 

19: 

Nielson  MlnNorm 

.0704 

.0100 

21: 

Hardy  Quadric 

.0397 

.00570 

23: 

Duchon  -  TPS 

.0588 

.00810 

27: 

Hardy  Reclp.  Quad. 

.0443 

.00528 

30: 

Foley  III 

.0823 

.00853 

Deviations  from  Saddle  test  surface,  25  points 


Table  D.3.3 


. 1  .".ii  j ; 
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RMS 

Deviation 

.0171 

.0202 

.0203 

.0189 

.0183 

.0189 

.0171 

.0205 

.0172 

.00952 

.0137 

.00955 

.0165 
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Method 

Maximum 

Deviation 

Mean 

Deviation 

1- 

Franke  -  3 

.0247 

.00491 

4 

Aklma 

.0256 

.00541 

10 

Aklma  Mod.  I 

.0248 

.00541 

13 

Nielson  -  Franke  Q 

.0340 

.00562 

14 

Mod.  Quad.  Shepard 

.0227 

.00529 

16 

Aklma  Mod.  Ill 

.0232 

.00575 

24 

Franke  -  TPS 

.0245 

.00440 

28 

Lawson 

.0234 

.00399 

19:  Nielson  MlnNorm 

.0161 

.00307 

21 

Hardy  Quadric 

.00709 

.00107 

23 

Duchon  -  TPS 

.0128 

.00265 

27 

Hardy  Reclp.  Quad. 

.00528 

.00055 

a  A  a  f* 

30:  Foley  III 

.0224 

%  00436 

Deviations 

from  Gentle  test  surface, 

25  points 

RMS 

Deviation 

.00651 

.00686 

.00681 

.00746 

.00669 

.00760 

.00556 

.00541 

.00433 

.00158 

.00351 

.00089 

.00588 


Method 


1 :  Franke  -  3 
4:  Aklma 
10:  Aklma  Mod.  I 
13:  Nielson  -  Franke  Q 
14:  Mod.  Quad.  Shepard 
16:  Aklma  Mod.  Ill 
24:  Franke  -  TPS 
28:  Lawson 

19:  Nielson  MlnNorm 
21:  Hardy  Quadric 
23:  Duchon  -  TPS 
27:  Hardy  Reclp.  Quad. 
30:  Foley  III 


Maximum 

Deviation 

.113 

.0534 

.0520 

.0550 

.0468 

.0510 

.0317 

.0455 

.0314 

.0189 

.0233 

.0144 

.0743 


Mean 

Deviation 

.0178 

.0108 

.0103 

.00890 

.00911 

.00908 

.00756 

.0129 

.00487 

.00453 

.00462 

.00288 

.0107 


RMS 

Deviation 

.0257 

.0149 

.0140 

.0127 

.0126 

.0128 

.0100 

.0286 

.00694 

.00595 

.00653 

.00386 

.0161 


Deviations  from  Steep  test  surface,  25  points 


Method 


1:  Franke  -  3 
4:  Aklma 
10:  Aklma  Mod.  I 
13:  Nielson  -  Franke  Q 
14:  Mod.  Quad.  Shepard 
16:  Aklma  Mod.  Ill 
24:  Franke  -  TPS 
28:  Lawson 

19:  Nielson  MlnNorm 
21:  Hardy  Quadric 
23:  Duchon  -  TPS 
27:  Hardy  Reclp.  Quad. 
30:  Foley  III 


Deviations 


Maximum 

Mean 

Deviation 

Deviation 

.0323 

.00498 

.0646 

.00903 

.0634 

.00811 

.0174 

.00199 

.0190 

.00200 

.0231 

.00190 

.0482 

.00690 

.0212 

.00216 

.0412 

.00470 

.0371 

.00403 

.0581 

.00557 

.0628 

.00774 

.0305 

.00568 

from  Sphere  test  surface,  25  points 
Table  D.3.6 


RMS 

Deviation 

.00748 

.0132 

.0121 

.00324 

.00336 

.00303 

.0106 

.00378 

.00765 

.00650 

.00925 

.0123 

.00822 


-111 


k 

V 


Test  Surface 

1 

2 

3 

4 

5 

6 

Data  Set 

100  points 

16 

28 

14 

16 

13-16-28 

13 

33  points 

13 

10 

14 

16 

13 

28 

25  points 

24 

10 

1-24 

28 

24 

16 

Local  Method 

With  Smallest  Deviation 

Table  E.l 

Test  Surface 

1 

2 

3 

4 

5 

6 

Data  Set 

100  points 

21 

21 

21 

21 

21 

13 

33  points 

27 

10 

21 

21 

21 

28 

25  points 

31 

30 

21 

21 

27 

16 

Method  With  Smallest  Deviation 


Method 

Parameter  Value 

Max. Dev. 

Mean  Dev. 

RMS  Dev 

1: 

Franke  -  3 

NPPR  ■  4 

.144 

.0113 

.0210 

6 

.0919 

.00842 

.0148 

8 

.0831 

.00864 

.0145 

V. 

Aklma 

NCP  -  4 

.130 

.00857 

.0153 

6 

.0647 

.00787 

.0125 

8 

.0934 

.00925 

.0156 

10: 

Akima  Mod,  I 

NCP  -  4 

.152 

.00898 

.0169 

6 

.0856 

.00784 

.0133 

8 

.0696 

.00874 

.0146 

13: 

Nielson  -  Franke  Q 

NPPR  -  12 

.0997 

.00729 

.0126 

18 

.0782 

.00741 

.0122 

24 

.0899 

.00831 

.0139 

14: 

Mod. Quad. Shepard 

NPPR  -  6  -  12 

.0663 

.00704 

.0117 

9  -  18 

.0573 

.00785 

.0128 

12  -  24 

.0735 

.00894 

.0148 

16: 

Aklma  Mod.  Ill 

NCP  ■  12 

.101 

.00709 

.0124 

18 

.0520 

.00729 

.0117 

24 

.0599 

.00821 

.01 33 

21: 

Hardy  Quadric 

NPPR  «  15 

.0287 

.00303 

.00578 

25 

.0225 

.00181 

.00357 

35 

.0185 

.00138 

.00257 

24: 

Franke  -  TPS 

NPPR  -  4 

.146 

.0104 

.0203 

6 

..0940 

.00887 

.0164 

8 

.0919 

.00804 

.0150 

27: 

Hardy  Reclp.  Quad. 

NPPR  »  15 

.0912 

.00601 

.0129 

25 

.0247 

.00283 

.00518 

35 

.0220 

.00217 

.00399 

30: 

Foley  III 

NIT  *  1 

.104 

.00745 

.0155 

3 

.0636 

.00473 

.00941 

5 

.0449 

.00376 

.00707 

Deviations  from  Exponential  test  surface,  100  points,  varying  parameters 


Method 

Parameter  Value 

Max. Dev. 

Mean  Dev. 

RMS  Dev 

1: 

Franke  -  3 

NPPR 

«  4 

.600 

.0446 

.0775 

6 

.240 

.0359 

.0486 

9 

.234 

.0428 

.0595 

4: 

Aklma 

NCP 

■  4 

.133 

.0256 

.0369 

6 

.134 

.0282 

.0386 

8 

.153 

.0302 

.0430 

10: 

Aklma  Mod.  I 

NCP 

■  4 

.133 

.0255 

.0371 

6 

.129 

.0280 

.0390 

8 

.146 

.0301 

.0432 

13: 

Nielson  -  Franke  Q 

NPPR 

■  12 

.214 

.0394 

.0670 

18 

.153 

.0350 

.0478 

24 

.132 

.0322 

.0433 

14: 

Mod.  Quad.  Shepard 

NPPR 

■  6  - 

12 

.230 

.0372 

.0549 

9  - 

18 

.158 

.0353 

.0486 

i 

12  - 

24 

.135 

.0338 

.0456 

16: 

Aklma  Mod.  Ill 

NCP 

-  12 

.176 

.0394 

.0560 

18 

.155 

.0355 

.0484 

24 

.127 

.0319 

.0433 

21: 

Hardy  Quadric 

NPPR 

■  15 

.120 

.0225 

.0307 

25 

.119 

.0235 

.0322 

35 

.129 

.0280 

.0397 

24: 

Franke  -  TPS 

NPPR 

-  4 

.186 

.0318 

.0455 

6 

.129 

.0267  ■ 

.0374 

9 

.143 

.0281 

.0404 

27: 

Hardy  Reclp.  Quad. 

NPPR 

■  15 

.122 

.0239 

.0333 

25 

.119 

.0214 

.0294 

35 

.119 

.0234 

.0323 

30: 

Foley  III 

NIT  « 

*  1 

.191 

.0238 

.0355 

3 

.165 

.0196 

.0310 

5 

.154 

.0209 

.0324 

Deviations  from  Exponential  test  surface,  25  points,  varying  parameters 

Table  F.3 
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Method 

Preprocessing 

Evaluation 

Total 

1: 

Frank®  -  3 

1.1 

8.8 

9.9 

4: 

Aklma 

2.2 

0.8 

3.0 

10: 

Aklma  Mod.  I 

2.8 

0.8 

3.6 

13: 

Nielson  -  Franke  Q 

10. 

1.9 

12. 

14: 

Mod. Quad.  Shepard 

8.6 

15. 

24. 

16: 

Aklma  Mod.  Ill 

11. 

0.8 

12. 

24: 

Franke  -  TPS 

2.7 

6.5 

9.2 

28: 

Lawson 

1.8 

1.7 

3.5 

19: 

Nielson  MlnNorm 

6.7 

3.8 

9.5 

21: 

Hardy  Quadric 

7.1 

13. 

20. 

23: 

Duchon»TPS 

7.5 

17.8 

24. 

27: 

Hardy  Reclp.  Quad. 

7.1 

13. 

20. 

30: 

Foley  III 

15. 

11. 

26. 

2 

Mod.  Shepard  8  Plane 

2.1 

25. 

27. 

3 

Mod.  Linear  Shepard 

1.2 

15. 

16. 

5 

McLain  M10 

— 

no. 

no. 

6 

Franke  -  1 

1.0 

8.0 

9.0 

7 

Mod.  Shepard 

— 

12. 

12. 

8 

Mod.  McLain  Mg 

— 

14. 

14. 

n 

Aklma  Mod.  II 

2.8 

.8 

3.6 

12 

Nielson  -  Franke  L. 

2.4 

1.5 

3.9 

17 

Quad.  Shepard 

33. 

22. 

55. 

18 

Shepard 

— 

17. 

17. 

20 

Rotated  Gausslans 

7.1 

13. 

20. 

22 

Duchon 

7.4 

15. 

22. 

25 

Foley  I 

— 

13. 

13. 

26 

Foley  II 

4.0 

16. 

20. 

29 

Rotated  B-Spllnes 

7.7 

23. 

31. 

31 

Foley  IV 

11. 

6.3 

17. 

Timing  :  100  points 


il 


■f'/' 

J 

1 

1 

1 

1  c 

Method 

Preprocessing 

Evaluation 

Total 

¥r' 

1: 

Franke  -  3 

0.3 

7.8 

8.1  1 

4: 

Aklma 

0.5 

0.7 

1.2  3 

q 

10: 

Aklma  Mod.  1 

0.7 

0.7 

1.4 

|i 

13: 

Nielson  -  Franke  Q 

2.4 

1.8 

4.2 

14: 

Mod.  Quad.  Shepard 

2.1 

5.6 

7.7 

1! 

16: 

Aklma  Mod.  Ill 

2.7 

0.7 

3.4 

ii 

24: 

Franke  -  TPS 

0.6 

4.6 

5.2 

1 

28: 

Lawson 

0.5 

1.5 

2,0 

i 

19: 

Nielson  MlnNorm 

1.9 

3.1 

5.0 

p 

21 : 

Hardy  Quadric 

0.5 

4.0 

4.5 

6't  ■ 

23: 

Duchon  TPS 

0.5 

5.3 

5.8 

27: 

Hardy  Reclp.  Quad. 

0.5 

4.0  ' 

4.5 

1  ; 

30: 

Foley  III 

1.6 

4.0 

5.6 

1 

2: 

Mod.  Shepard  9  Plane 

0.2 

9.9 

10. 

3: 

Mod.  Linear  Shepard 

0.2 

5.1 

5.3 

ti  ■ 

5: 

McLain  M^ 

— 

50. 

50. 

6: 

Franke  -  1 

0.3 

6.9 

7.2 

7: 

'Mod.  Shepard 

■  ■  • 

4.6 

4.6 

1?  ( 

8: 

■Mod.  McLain  Mg 

— 

5.7 

5.7 

11: 

Aklma  Mod.  II 

0.7 

0.7 

1.4 

12: 

Nielson  -  Franke  L. 

0.3 

1.5 

1.8 

17: 

Quad.  Shepard 

4.1 

7.1 

11. 

18: 

Shepard 

6.4 

6.4 

20: 

Rotated  Gausslans 

0.5 

4.0 

4.5 

;■ 

22: 

Duchon 

0.5 

5.0 

5.5  i 

25: 

Foley  I 

... 

4.0 

4.0  [ 

26: 

Foley  II 

0.9 

8.2 

9.1  j 

29: 

Rotated  B-Splines 

0.5 

7.8 

8.3 

i 

31 : 

Foley  IV 

1.1 

2.7 

3.8  j 

Timing: 

33  points 

Table  T. 2 

C 

I 

-118- 

|GW.. 

'  i  HuU-'1  « Ilk Hjiv ieJtt ; -  tv;  fV'4'r 

Method 

1:  Franke  -  3 
4:  Akima 
10i  Akima  Mod.  1 
13:  Nielson  -  Franke  Q 
14:  Mod.  Quad.  Shepard 
16:  Akima  Mod.  Ill 
24:  Franke  -  TPS 
28:  Lawson 

19:  Nielson  MlnNorm 
21:  Hardy  quadric 
23:  Duchon  -  TPS 
27:  Hardy  Reclp.  Quad. 
30:  Foley  III 


Preprocessing 

0.2 

0.3 

0.5 

1.7 
1.5 

1.8 
0.4 
0.4 

1.4 

0.2 

0.2 

0.2 

1.1 


Evaluation  Total 


7.7 
0.6 
0.6 

1.7 
4.2 
0.6 
4.5 
1.4 


7.9 
0.9 
1.1 

3.4 

5.7 

2.4 

4.9 

1.8 


2.9 

3.1 

4.0 

3.1 

3.1 


4.3 

3.3 

4.2 

3.3 
4.2 


2:  Mod.  Shepard  ®  Plane 
3:  Mod.  Linear  Shepard 
5:  McLain  M^g 

6:  Franke  -  1 
7:  Mod.  Shepard 
8:  Mod.  McLain  Ma 

11:  Akima  Mod.  II 

12:  Nielson  -  Franke  L. 

17:  Quad.  Shepard 

18:  Shepard 

20:  Rotated  Gausslans 

22:  Duchon 

25:  Foley  I 

26:  Foley  II 

29:  Rotated  B-Spllnes 

31:  Foley  IV 


0.1 

0,1 

0.2 


0.5 

0.2 

2.5 


0.2 

0.2 


0.7 

0.3 

0.8 


7.6 
4.0 

40. 

6.7 

3.7 

4.4 

0.6 

1.4 

5.5 

4.1 

3.1 

3.8 
3.0 
7.3 

5.9 

1.8 


7.7 

4.1 
40. 

6.9 

3.7 

4.4 

1.1 
1.6 
8.0 

4.1 
3.3 
4.0 
3.0 
8.0 

6.2 
2.6 


Timing  :  25  Points 

Table  T.3 
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Foley's  Iterated  Modified  Shepard 
Delta  Sum  Bicubic  Spline 

Figure  1.4.1.31 
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Franke's  Method,  Mode  ■  3 
Figure  1 .4.4.1 
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Hardy  Reciprocal  Multiquadric  Method 
Figure  1 .4.4.27 
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Nlel son-Franke  Quadratic  Triangular  Method 
Figure  1.4.5.13 
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Nlelson-Franke  Quadratic  Triangular  Method 
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Rotated  Cubic  B-Spllne  Basis 
Figure  K.O.o.29 
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Figure  3.3.1 .16 


0  PT  CLIFF  FUNCTION  30378  1054 

TEST  SURFACE  PLOT 

Figure  4, 0.2.0. 


244 


5%-Sr-' 

5gSSsa5r 


0  PT  SADDLE  MOUNTAIN  THREE  31778 

113 

1 

TEST  SURFACE  PLOT 

Figure  4. 0.3.0. 

•245- 

■  '  . . . . . .  .  •-•  - . -  . . . .  .........  ,  ..  _ _ _ _ _ _ _ 

ANW^WWiittsiuii: 

100PT  EXPONENTIAL  HUMPS,  DIP 


C  HARDY  MQ  SURFACE 


NPPR  = 
Figure  4.1 .1 .21 


1 GOPT  EXPONENTIAL.  HUMPS.  DIP  8087S 
FRANKE  W/  THIN  PLATE  NPPR  =  6 


Figure  4.1.1.24 


100PT  EXPONENTIAL  HUMPS,  DIP  41079  1825 

HARDY  RECIP  MQ  SURFACE  NPPR  =  25  B  =  0.185 


Figure  4.1.1.27 


Figure  4.1 .1 .28 


100PT  EXPONENTIAL  HUMPS,  DIP  S0P79 


FOLEY  ITERATED  B I  CUB  ICS  NPPR 

Figure  4.1 .1.30 


100  PT  CLIFF  FUNCTION  30878  1121 

'RANKE  -  SARD  BS  PLANE  NPPR  =  S 


Figure  4.1 .2.1 


I  100PT  CLIFF  FUNCTION 

t  HflRDT  MQ  SURFACE  NPPR  = 

(Figure  4.1 .2.21 

C 

■  1^  ■■l-Vh fi'iM-1 '  ■?’>! t »•> .1  mUifcjftfilirtwn  ^ ...  I "  .i.n.  I 


31579  2015 

25  R  =  0.185 


LW'i'USil.. 


1O0PT  CLIFF  FUNCTION  40679  1525 

HARDY  RECIP  MQ  SURFACE  NPPR  =  25  R  =  0.165 


Figure  4.1 .2.27 


100  PT  SADDLE  MOUNTAIN  THREE  50278  1502 

FRflNKE  -  SARD  BS  PLANE  NPPR  =  6 


Figure  4,1.3.10 


.  . . .  .  »«M~. . ....  ■> 


100PT  SADDLE  MOUNTAIN  THREE 


> 


1 


DUCHON  S  THIN  PLATE 

c 


Mi*  t  i"l  it;.5!  «i»t wr 1 1  •- 


■ulitfv 


Figure  4.1.3.23 


-.-.f  -.I.H.  .-I.I 


1Q0PT  SADDLE  MOUNTAIN  THREE  80373 

FRANKE  W/  THIN  PLATE  NPPR  --  6 

Figure  4.1.3.24 


SMnU'U.II  If.'fil'I  l"  ill!.-  ■ 


100PT  SADDLE  MOUNTAIN  THREE  ,  80279  i  1 01 

FOLEY  ITERATED  BICUBICS  NPPR  =  3 

Figure  4.1.3.30 

.  ......  i  .1  .  -I  .Ifll-iiVt..,  -..itt.a..  i.-.J.  .  !.  f..,..  -.1  ..  11! ...  .  .i  .wi-.  •-  ^  <•••  I  -.ft  ■--.■■■■.-..-J-.-.  WIM-I-Oi.  I. 


z 


33PT 

EXP 

HUMPS. 

DIP  ON  S P hI F i 3 E 

4  1  8  79 

1  02  1 

RKlMfl  - 

FEM 

F  I  T 

Nif-'PR 

Figure  4.2.1 .4 

P 

REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


33PT  EXPONENTIAL  ON. SPARSE  1  1979  1  208 

NIELSON-FRANKE  QUAD-TRI  NPPR  -  18  R  =  0.522 


Figure  4.Z.1 .13 


33PT  EXPONENTIAL  ON  SPARSE  11979  1202 

QUADRATIC  SHEPARD  METHOD  NPPR  =  918  R  =  0.369 


Figure  4.2.1.14 


33PT 

NIELSON 


EXP  HUMPS,  DIP  ON  SPRRSE  12979  1351 

MIN  NORM  NETWORK 


Figure  4.2.1.19 


33PT  EXP  HUMPS.  PIP  ON  SPARSE  80379  910  * 

i 

FRANKE  W/  THIN  PLATE  NPPR  =  6 

(  Figure  4.2.1.24 

REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


33PT  EXP  HUMPS,  DIP  ON  SPRRSE 


HARDY  RECIP  MQ  SURFACE 


NPPR  = 


Figure  4.2.1 .27 


3 3 P T  CLIFF  FUNCTION  ON  SPARSE 
FRfiNKE  -  SARD  BS  PLANE  NPPR  = 

Figure  4. 2. 2.1 


•*+*#$**  ** 


KfZSgnA 


MijMl 

7Mrt  1}  ini 


lmjj$  J  j 

^<S^yinjjifUkriU  \  if 

^Q^ynilUrfiii  a,]  :  r 

"- t^s.7  -.-*  ”  /  /  /  .j  .'  ,  i  > 


^irr  } 

NJ.  ll- 


33PT  CLIFF  FUNCTION  ON  SPHF1SE  'iiO/U 


IMfl  -  FEM.  F  IT 


NPPFi 

Figure  4. 2. 2. 4 


REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


..uWlsiiii-iJtfci-.  ■i-»-  -«■  ^^'-^liiftiiiM^l  fKi'lin  II  -  CjTi  __ ..  tj.  . 


3 3 P T  CLIFF  FUNCTION  ON  SPARSE  31579 
RDY  MQ  SURFACE  NPPR  -  25  R  = 


2017 

0.308 


Figure  4.2.2.21 


33PT  CL.  I FF  FUNCTION  ON  SPARSE 


803  79 


3  3PT 
AKIMA  - 


SADDLE  ON  SPARSE  41879 

FEM  FIT  NPPR  -  6 

Figure  4. 2. 3. 4 


REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


I 


33PT 

SADDLE 

ON  SPARSE 

31079 

1909 

NIELSON 

-FRANKE 

QURD-TR I  NPPR  = 

18  R  = 

0.522 

Figure  4. 2. 3.1: 

3 

..  •  — . . j 

. . . 

. . -fc- . . . . 

33PT  SADDLE  ON  SPARSE 


31179 


1441 


QUADRATIC  SHEPARD  METHOD  NPPR  =  918  R  =  0.369 


Figure  4.2.3.14 


DUCHCN  S  THIN  PLRTE 


Figure  4.2.3.23 


33PT  SADDLE  ON  SPARSE 


30879 


FRANKE 

c 


W/  THIN  PLATE  NPPR  = 

Figure  4.2.3.24 


REPRODUCED  FROM 

BEST  AVAILABLE  COPY 


3 3 P T  SADDLE  ON  SPARSE  32379 

LAWSON  TRIANGLE  METHOD 

Figure  4.2.3.28 


1422 


L  1.1. f I V.  iH:  V..1  l.-J'Jl  1.1 -11  M.lit.'IM.' IV  ■ .  t  ITU;  ,1  .1-  .  j|lfl  •'fe.-l,  (hjftdj.'filfctf  I'JAiis  .'i  i  t&i  A.fttUu  < 


33PT  SADDLE  ON  SPARSE 


FOLEY  ITERATED  BICUBJCS  NPPR 

Figure  4.2.3.30 


DIP  ON  GREGS  4  I  ft  7  9 

NPF’R  =  G 

Figure  4.3.1 .4 

REPRODUCED  FROM  -325 

BEST  AVAILABLE  COPY 


25PT  EXPONENTIAL  ON  GREGS  11579  1209 

QUADRATIC  SHEPARD  METHOD  NPPR  =  %tB  R  =  0.414 

Figure  4.3.1 .14 

. _  ......  .  ... 


r*  aw  «.■„  'ini'  w  *  *  -  u-sv .,  ..<i.*yfe»i.- *.•••• . ^h-m^  4wMu,w'r*v/ 


I WW *&H n if  tarrr,  fa^va  imvJLUmm *fcl  ft»W» ~<.VW-rv  » 


z 


?5PT  EXP  HUMPS,  DIP  ON  GREGS  20679  1535 

DUCHON  S  THIN  PLATE 


Figure  4.3.1.23 


25PT  EXP  HUMPS,  DIP  ON  GREGS 
LAWSON  TRIANGLE  METHOD 


Figure  4.3.1 .28 


25PT  EXP  HUMPS,  DIP  ON  GREGS 
FOLET  ITERATED  BICUBICS  NPPR  = 

Figure  4.3.1.30 


60279 


3 


REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


25PT  CLIFF  ON  GREGS 
fiKIMfl  -  FEM  FIT  NPPR  = 

Figure  4. 3. 2. 4 


41879  1022 

6 


REPRODUCED  FROM 
BEST  AVAILABLE  COPY 


REPRODUCED  FROM 


25PT 

SRDDLE 

ON  GREGS 

31079 

NIELSON 

-FRflNKE 

QUflD-TRI  NPPR  = 

18  R  = 

a'i,y.|i.iiah'au>''’',''iirH".u'. 

Figure  4.3.3.13 

25PT  SADDLE  ON  GREGS 


80879 


FRflNKE  W/  THIN  PLATE  NPPR  = 

Figure  4,3.3.24 


6 


Iri'.i  tUUftisiill Jittv i*  i  MhiAfi 1 1  itffflilfct Miti'id ''iJtfiiAil Kr i' (  ■iii  1  il Jillidlt IM '< V UiUB i i ?lf j  'inti; 


tofctilrcflSWWA'  '  '"'hTr'^iVii-  -tfiiAilaiir uu 


►***.  '1 


— n 


DISTRIBUTION  LIST 


Defense  Documentation  Canter  2 

Cameron  Station 
Alexandria,  VA  22314 

Dudley  Knox  Library  2 

Navel  Postgraduate  School 
Monterey*  CA  93940 

Dean  of  Research  2 

Naval  Postgraduate  School 
Monterey*  CA  93940 

Department  of  Mathematics 

C.  0.  Wilde*  Chairman  1 

F.  D.  Faulkner,  Acting  Chairman  1 

A.  L.  Schoenstadt  1 

R.  Franke  15 

Naval  Postgraduate  School 
Monterey,  CA  93940 

Dr.  Richard  Lau  1 

Office  of  Naval  Research 
'030  East  Green  St. 

Pasadena,  CA  91106 

Professor  R.  E.  Barnhill  1 

Department  of  Mathematics 
University  of  Utah 
Salt  Lake  City.  UT  84112 

Professor  G.  M.  Nielson  1 

Department  of  Mathematics 
Arizona  State  University 
Tempe,  AZ  85281 

Chief  of  Naval  Research  2 

ATT:  Mathematics  Program 
Arlington,  VA  22217 

Rosemary  E.  Chang  1 

Sandla  Laboratories 

Applied  Mathematics.  Division  2-8235 

Livermore,  CA  94550 


Hiroshi  Aklma 

Office  of  Telecommunications 
Department  of  Commerce 
Boulder,  CO  80302 


1 


MISSING 

IN 

ORIGINAL 

DOCUMENT 


BPNWPMWIWW  B**^***!^^  ..rw-  rt******  «c*  ■ 


Eric  A.  Lunds trom 
Cod*  3275 
Navtl  Weapons  Center 
China  Lake.  CA  93555 

John  R.  Nice 
Math  Science  428 
Purdue  University 
W.  Lafayette,  IN  47907 

Or.  Robert  Rlffenburgh 
Code  18 

Naval  Ocean  System  Center 
San  Diego,  CA  92152 

John  K.  Munro,  Jr. 

Bldg.  9104-2 

Oak  Ridge  National  Laboratory 

P.  0.  Box  Y 

Oak  Ridge,  TN  37830 

Paul  Peterson 
Woodward  Governor  Co. 

5001  N.  2nd  St. 

Rockford.  I L  61101 

John  Bobbitt 
Continental  Oil  Co. 

P.  0.  Box  1267 
Ponca  City,  OK  74601 

Eric  Grosse 

Computer  Science  Department 
Stanford  University 
Stanford,  CA  94305 

LCDR  Alan  J.  Plckrell 
3233  NE  105th  Street 
Seattle,  WA  98125 

Dr.  Farhad  Rajabl 
SMC  #1717 

Naval  Postgraduate  School 
Monterey,  CA  93940 

Visiting  Professor  0.  C.  Zlenklewlcz 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  CA  93940 

Arthur  R.  Paradis 
President,  Dynamic  Graphics,  Inc. 
2150  Shettuck  Avenue 
Berkeley,  CA  94704 


Or.  Paul  T.  Boggs 
Army  Research  Office 
P.  0.  Box  1211 
Triangle  Park,  NC  27709 

Stavros  Busenberg 
Department  of  Mathematics 
Harvey  Mudd  College 
Claremont,  CA  91711 

Or.  Theodore  Rubin 

IBM  Watson  Research  Center 

P.  0.  Box  218 

Yorktown  Heights,  NY  10598 

Professor  J.  Ben  Rosen 
Computer  Science  Department 
University  of  Minnesota 
Minneapolis,  MN  5*5455 

Andy  Schoene 

Computer  Science  Department 
General  Motors  Research  Laboratory 
Warren,  MI  48090 

Professor  L.  L.  Schumaker 
Department  of  Mathematics 
University  of  Texas  at  Austin 
Austin.  TX  78712 

L.  F.  Shamplne 

Numerical  Mathematics  Dlv.  5122 
Sandla  Laboratories 
Albuquerque,  NM  87185 

Patrick  W.  Gaffney 
Union  Carbide/Nuclear  Division 
Computer  Sciences  Division 
Box  Y 

Oak  Ridge,  TN  37830 

Dr.  Carl  de  Boor 
Mathematics  Research  Center 
University  of  Wisconsin 
Madison,  WI  53706 

Grant  Burgt's 
Mobil  Oil  Co. 

1001  Howard  Avenue,  Room  937 
New  Orleans,  LA  70113 


Jean  Duchon 
University  of  Grenoble 
F- 38041 

Grenoble.  FRANCE 

J.  K.  Reid 
AERE 

Building  8-9 

Harwell.  Oldcot.  Blrkshlre.  ENGLAND 

Professor  John  A.  Soulier 
Department  of  Mathematics 
North  Carolina  State  University 
Raleigh.  NC  27607 

G.  Bolondl 
AGIP 

S.  Donato  Milanese 
Milan.  ITALY 

D.  H.  McLain 
Computing  Centre 
University  of  Sheffield 
Sheffield.  United  Kingdom  S10  2TN 

Nell  Stahl 

U.  w.  -  Fox  Valley 

Midway  Rd. 

Menuyha,  HI  54952 

C.  A.  Steele,  Jr. 

P.  0.  Box  45 
Magnolia.  MA  01930 

Kunlo  Tanabe 
Applied  Math  Department 
Brookhaven  National  Laboratory 
Upton,  NY  11973 

Mr.  B.  Clifford 

Department  of  Computer  Services 
University  of  Calgary 
Calgary,  Alberta,  CANADA  T2N  1N4 

James  R.  Jancaltus 
USA  Topographic  Labs 
Ft.  Bel voir,  VA  22060 

John  l.  Junklns 

Department  of  Engineering  Science 

Virginia  Polytechnic  Institute  and  State  University 

Blacksburg,  VA  24061 


Harry  Joseph  Feeney  III 
Naval  Ocean  Systems  Center 
San  01 ego.  CA  92152 

Dave  Cooper 
AFHRL/ASM 

Wrlght-Patterson  AFB,  OH  45433 

Stanley  C.  Elsenstat 
Department  of  Computer  Science 
Yale  University 
10  Hlllhouse  Avenue 
New  Haven,  CT  06520 

Professor  Grace  Wahba 
Department  of  Statistics 
University  of  Wisconsin 
Madison,  HI  53706 

Jesse  Y.  Wang 

Applied  Mathematics  Division 
Argonne  National  Laboratory 
9700  So.  Cuss  Avenue 
Argonne,  IL  60439 

Shraya  Yosefa 

Department  of  Applied  Mathematics  and  Statistics 
State  University  of  New  York  at  Stony  Brook 
Stony  Brook,  NY  11794 

R.  E.  Funderllc 

P.  0.  Box  X,  Bldg.  4500N,  0224 
Oak  Ridge,  TN  37830 

A.  J.  Goldman 

Chief,  Operations  Research  Dlv. 

U.  S.  Department  of  Commerce 
National  Bureau  of  Standards 
Washington,  D.  C.  20234 

Dr.  Gary  Herron 
Boeing  Computer  Services 
P.  0.  Box  24346 
Seattle,  WA  98124 

Dr.  R.  Peter  Dube 
Boeing  Computer  Services 
P.  0.  Box  24346 
Seattle,  WA  98124 

Professor  Gerald  D.  Taylor 
Department  of  Mathematics 
Colorado  State  University 
Ft.  Collins,  CO  80523 


1 


Professor  Tow  Foley 
Department  of  Computer  Science 
California  Polytechnic  State  University 
San  Luis  Obispo,  CA  93407 

William  L.  Vlttltow 
Technical  Services 
Lockheed  -  California 
P.  0.  Box  551 
Burbank,  CA  91503 

M.  J.  0.  Powell 

Department  of  Applied  Mathematics 
University  of  Cambridge 
Silver  Street 
Canfcrldge,  ENGLAND 

Jean  Melnguet 

Istltute  de  Mathematique  P.  et  A 
Unlverslte  de  Louvain 
Chemin  du  Cyclatron  2 
B-1348  Louvain- la-Nerlve 
BELGIUM 

A.  R.  Forrest 

School  of  Computing  Studies 
University  of  East  Anglia 
Norwich,  ENGLAND  NR4  7TJ 

Professor  Gordon  E.  Latta 
Department  of  Mathematics 
Naval  Postgraduate  School 
Monterey,  CA  93940 

Professor  Garrett  Blrkhoff 
Mathematics  Department 
One  Oxford  Street 
Cambridge,  MA  02138 

Commanding  Officer 

Fleet  Numerical  Oceanographic  Center 
Naval  Postgraduate  School 
Monterey,  CA  93940 


Officer  in  Charge  _  .... 

Naval  Environmental  Research  Prediction  Facility 
Naval  Postgraduate  School 
Monterey,  CA  93940 


j  q  |  H3V6S 

Division  of  Numerical  Analysis  and  Computing 
National  Physical  Laboratory 
Teddington,  ENGLAND  TW11  OLW 


i 


i 


i 


i 


i 


l 


36? 

•i 


Dr.  David  Kahaner 
U.  S.  Department  of  Commerce 
National  Bureau  of  Standards 
Scientific  Computing  Division 
Washington,  D.  C.  20234 

Donald  Beardsley 
Retro  Lewis  Corporation 
One  Energy  Center 
P.  0.  Box  2250 
Denver,  CO  80201 

Dr.  James  G.  Smith 
Office  of  Naval  Research 
Arlington,  VA  22217 

Frank  Hagln 
Computing  Services 
University  of  Denver 
Denver,  CO  80208 

Professor  Roll  and  L.  Hardy 
Department  of  Civil  Engineering 
Iowa  State  University 
Ames,  IA  50011 

John  M.  Karon 

Department  of  Biostatistics 
School  of  Public  Health 
Trailer  #32  306H 
University  of  North  Carolina 
Chapel  Hill,  NC  27514 

R.  W.  Klopfensteln 
RCA  Laboratories 
Princeton,  NJ  08540 

Gene  Golub 

Computer  Science  Department 
Stanford  University 
Standord,  CA  94305 

Professor  William  J.  Gordon 
Department  of  Mathematics 
Drexel  University 
Philadelphia,  PA  19104 

Charles  Mlcchelll 

Thomas  J.  Watson  Research  Center 

Post  Office  Box  218 

Yorktown  Heights,  NY  10598 


Michael  Mlnkoff 
Argonne  National  Laboratory 
9700  South  Cass  Avenue 
Argonne.  IL  60439 

R.  K.  Rew 
NCAR 

P.  0.  Box  3000 
Boulder.  CO  80307 

H.  Flchtner 
Bell  Labs 

Murray  Hill.  NJ  07974 

Peter  R.  Eiseman 
I  CASE 
MS  132C 

NASA  Langley  Research  Center 
Hampton,  VA 

Stephen  C.  Banks 
Sun  Production  Co. 

503  N.  Central  Expressway 
Richardson,  TX  75080 

Pablo  Barrera 
Facultao  de  Cl encl as 
Department  de  Matematlcas 
U.  N.  A.  M. 

Mexico  20.  0.  F.  MEXICO 

Richard  A.  Hansen 
308  TMCB 

Brigham  Young  University 
Provo,  UT 

P.  H.  Merz 
Chevron  Research  Co. 
p.  0.  Box  1627 
Richmond,  CA  94802 

C.  Bunch 
P.  0.  Box  1267 
Conoco  R  &  0 
Ponca  City,  OK  74601 

James  Lyness 
AMD-ANL 

Argonne,  II  60439 

Fred  N.  Fritsch 
Lawrence  Livermore  Laboratory 
P.  0.  Box  808  (L-300) 
Livermore,  CA  94500 


W.  J.  Schaffers 

Dupont  Co.,  Engineering  Department 
Exp.  Station,  304 
Wilmington,  DE  19898 

U.  Ascher 

Department  of  Computer  Science 
University  of  British  Columbia 
Vancouver,  BC  CANADA 

Joe  McGrath 
KMS  Fusion,  Inc. 

3941  Research  Park  Dr. 

Ann  Arbor,  MI  48104 

P.  S.  Jensen 

Lockheed  Research  5233/205 
3251  Hanover  St. 

Palo  Alto,  CA  94304 

Myron  Ginsberg 
Computer  Science  Department 
General  Motors  Research  Laboratory 
Warren,  MI  48090 

L.  Kratz 

Mathematics  Department 
Idaho  State  University 
Pocatello.  ID  83209 

Alan  Pierce 
Amoco  Production  Co. 

P.  0.  Box  591 
Tulsa,  OK  74102 

A.  K.  Cline 

Department  of  Computer  Science 
University  of  Texas  at  Austin 
Austin,  TX  78712 

W.  Roy  Wessel 
CDC 

7995  E.  Prentice  Ave. 

Englewood,  CO  80111 

J.  W.  Chalmers 
HA0/NCAR 
Box  3000 

Boulder,  CO  80307 

Richard  B.  Evans 
Ocean  Data  Systems,  Inc. 

6000  Executive  Blvd. 

Rockville,  MD  20852 


1 


1 


1 


1 


1 


1 


1 


1 


1 

( 


372 


1 


L.  H.  Seltelman 
Pratt  and  Whitney  Aircraft 
400  Main  St. 

East  Hartford,  CT  06108 

John  A.  Carpenter  1 

Oak  Ridge  National  Laboratories 

Bldg.  4500-N 

Room  E  208 

P.  0.  Box  X 

Oak  Ridge,  TN  37830 

James  C.  Ferguson  1 

Los  Alamos  Scientific  Laboratory 
Los  Alamos,  NM  87544 

Professor  David  Salinas  1 

Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  CA  93940 

Dr.  Nlra  Rlchter-Dyn  1 

Department  of  Mathematical  Sciences 
Tel -Aviv  University 
Tel -Aviv,  ISRAEL 


sssssss, 


